Estimating Group Differences in Network Models

using Moderation Analysis

Jonas Haslbeck
Psychological Methods Group, University of Amsterdam

Abstract
Statistical network models such as the Gaussian Graphical Model and the Ising model have
become popular tools to analyze multivariate psychological data sets. In many applications
the goal is to compare such network models across groups. In this paper I introduce a method
to estimate differences in network models across groups that is based on moderation analysis.
This method is attractive because it allows to make comparisons across more than two groups
within a single model, and because it is implemented for all commonly used cross-sectional
network models. Next to introducing the method, I evaluate the performance of the proposed
method and existing approaches in a simulation study. Finally, I provide a fully reproducible
tutorial on how to use the moderation method to compare a network model across three groups
using the R-package mgm.

1 Introduction
Statistical network models such as the Gaussian Graphical Model (GGM), the Ising model or Mixed
Graphical Models (MGMs) have become popular tools to analyze multivariate cross-sectional
datasets (Epskamp, Maris, Waldorp, & Borsboom, 2016; Epskamp, Waldorp, Mõttus, & Bors-
boom, 2018; Haslbeck & Waldorp, 2020; van Borkulo et al., 2014; Williams & Mulder, 2019). In
many of these applications one is interested in comparing such network models between two or
more groups. For example, Fritz, Fried, Goodyer, Wilkinson, and van Harmelen (2018) compared
the relations between resilience factors in a network model for adolescents who did or did not
experience childhood adversity; van Loo et al. (2018) compared the relations between depression
symptoms across various environmental and genetic risk factors; and van Borkulo et al. (2015)
compared the relations between depression symptoms before and after treatment.
There are already several methods available to estimate group differences in network models.
The Network Comparison Test (NCT; van Borkulo et al., 2017), which is currently implemented
for the GGM and the Ising model, uses a permutation test to compare pairs of groups. The Fused
Graphical Lasso (FGL; Danaher, Wang, & Witten, 2014), which is currently available only for
the GGM, compares two or more groups by applying a penalty term to group differences and
performing model selection across different penalties. Epskamp, Isvoranu, and Cheung (2020)
proposed a SEM-based method that uses iterative model search and pruning. Finally, one can
compare GGMs across two or more groups in a Bayesian framework, which make use of the Bayes
factor or threshold the posterior of edge-differences (Williams, Rast, Pericchi, & Mulder, 2019).
In the present paper I introduce an additional method to estimate group differences in network
models that is based on moderation analysis. Specifically, the grouping variable is included as a
categorical moderator variable, and group differences are determined by estimating the moderation
effects. This method is attractive because it allows to make comparisons across more than two
groups within a single model, and because it is implemented for all commonly used cross-sectional
network models. Next to introducing this method, I compare the performance of all discussed
methods in a simulation study, a comparison that is currently missing in the literature. Finally, I
provide a fully reproducible tutorial on how to estimate group differences in network models with
the moderation method using the R-package mgm (Haslbeck & Waldorp, 2020).
Contact: jonashaslbeck@[Link] | [Link]

1
2 Detecting Group Differences in Network Models
We first review existing approaches to estimate differences in parameters of network models across
groups and then introduce the moderation approach.

2.1 Existing Approaches

A popular method is the Network Comparison Test (NCT; van Borkulo et al., 2017), which is
based on a three-step procedure. First, the network model of choice is estimated separately on
both groups A, B and the difference between the parameter estimates in both groups of a parameter
serves as the test statistic. Second, for each parameter difference, a sampling distribution under
the null hypothesis of no difference is created by repeatedly assigning cases randomly to two groups
with the same sizes as A, B, and computing the difference in the parameter estimates. This yields
a sampling distribution that is centered at zero. Finally, using a specified significance value α, it
is determined whether the test statistic is significantly different from zero. An advantage of this
generic method is that it can be applied to essentially any model.
The Fused Graphical Lasso (FGL; Danaher et al., 2014) allows one to simultaneously estimate
GGMs in two or more groups. This is achieved by adding an additional `1 -penalty term for differ-
ences between groups to the already existing `1 -penalty of the graphical lasso (Friedman, Hastie,
& Tibshirani, 2008). This leads to small group differences being set to zero, while large enough
differences are included in the model. The R-package EstimateGroupNetwork (Costantini & Ep-
skamp, 2017) implements a version of the FGL in which the regularization parameters are selected
either with cross-validation or the Extended Bayesian Information Criterion (EBIC) (Foygel &
Drton, 2010). Currently, this method is only available to jointly estimate multivariate Gaussian
distributions, however, the general method could be extended to other models as well.
Epskamp et al. (2020) recently suggested a SEM-based approach to estimating group differences
that consists of the following three steps. First, one estimates a model in each separate data set
and removes non-significant partial correlations. Second, a pooled model is estimated in which
each parameter is included that was included in at least one of the individual models from the first
step. Third, equality constraints are freed in a step-wise fashion until the BIC does not improve
any further. The model at the end of this procedure is taken as the final model. This method is
available in the R-package psychonetrics (Epskamp, 2020).
Williams et al. (2019) proposed two Bayesian methods to test for differences in parameters
across groups. The first method uses a Bayes factor to compare the hypothesis that a given partial
correlation is the same in all groups vs. not. The second method computes the posterior of the
difference between partial correlations in two groups, and uses a threshold α on the posterior to
decide whether the difference is reliably different from zero. Both methods are implemented in the
R-package BGGM (Williams & Mulder, 2019). Recently these methods have been extended to
data consisting of continuous, ordinal and binary variables using a semi-parametric copula model.
Finally, to test whether partial correlations in GGMs are different across two groups one can
use traditional Frequentist methods. For example, one can obtain unbiased estimates of partial
correlations in both groups, and transform them with Fishers Z-transformation to obtain param-
eters that are approximately normally distributed (Fisher, 1915). Subsequently one can use the
Z-score of the difference together with a threshold α to test the null hypothesis that the difference
between the partial correlations is zero.

2.2 Moderation Approach

In this section, I show how to estimate group differences using moderation with categorical mod-
erators. I first illustrate the idea with a simple regression example, and then extend it network
models. Consider the simple regression problem:

Y “ βY X X ` ε, (1)

involving the variables X and Y , a parameter βY X and a Gaussian error term ε.

We now introduce an additional categorical variable G which indicates group membership and
takes on the values 1, 2 and 3. We are interested in whether the parameter βY X differs across the
three groups, and therefore add the interaction between variables X and G to the model:

2
 Y “ βY X X ` βY G2 X IpG “ 2qX ` βY G3 X IpG “ 3qX ` ε, (2)

where IpG “ kq is the indicator function for the group at hand being group k. For example, if G is
equal to 3, then the terms IpG “ 1q and IpG “ 2q are equal to zero, and IpG “ 3q is equal to one.
How do we obtain the main effect of X on Y in each group from the parameters in Equation
2? In Group 1, the main effect is equal to βY X ; in Group 2, it is equal to βY X ` βY G2 X ; and
in Group 3, it is equal to βY X ` βY G3 X . We see that by adding the grouping variable G as an
interaction term we obtained the main effect of X on Y for each group, and therefore also the
differences across groups. In this example the grouping variable G had 3 categories (groups), but
the grouping variable can also be binary or have more than 3 categories. Note that these main
effects relate to pairwise interactions between variables network models. For example, for the GGM
the main effects can be standardized to obtain partial correlations.
We can use the same principle to estimate group differences in network models. This is because
network models can be estimated with the pseudo-likelihood approach, in which one estimates the
conditional distribution of each node, and then combines the estimates obtained from the condi-
tional distributions to the full network model (Meinshausen, Bühlmann, et al., 2006). In the case
of the multivariate Gaussian distribution and its implied Gaussian Graphical Model (GGM) with
p variables, estimating the conditional distribution of a given variable Xi amounts to estimating p
multiple regression models:

Xi “ βXi X1 X1 ` ¨ ¨ ¨ ` βXi Xp Xp ` ε, (3)

in which each variable Xi , i P t1, 2, . . . , pu is predicted by all variables except itself.

Similarly to Equation 2 in the simple linear regression problem, we can add a grouping variable.
Here we choose a binary grouping variable G:

Xi “ βXi X1 X1 ` βXi G1 X1 IpG “ 1qX1 ` ¨ ¨ ¨ ` βXi Xp Xp ` βXi G1 Xp IpG “ 1qXp ` ε. (4)

As mentioned above, in the pseudo-likelihood approach each variable is predicted by all other
variables in a multiple regression model. This leads to two estimates for each main effect: for
example, we obtain the estimate βX1 X4 from the regression on X1 and the estimate βX4 X1 from
the regression on X4 . These two estimates are estimates of the same pairwise interaction between
variables X1 and X4 . These estimates can be standardized to obtain the partial correlation ρ1,4
in the p-variate multivariate Gaussian distribution (Epskamp et al., 2018). To arrive at a single
estimate we aggregate them using the AND-rule (take average) or the OR-rule (take average if
both estimates are nonzero, otherwise set to zero). Similarly, we obtain two estimates for the
interaction parameters which captures the group differences.
The same procedure can be applied to the Ising model. The only difference is that the nodewise
regressions are not linear but logistic regressions (Epskamp et al., 2016). More generally, the
procedure can be applied to Mixed Graphical Models (MGMs) (Yang, Baker, Ravikumar, Allen, &
Liu, 2014), which generalize both the GGM and the Ising model. In this setting, each variable is a
conditional exponential family distribution (e.g., Gaussian, Binomial, Poisson), and one can apply
the above procedure by estimating the appropriate regressions in the Generalized Liner Model
(GLM) framework. These models are generalized by higher-order (or moderated) MGMs, which
are described in detail in Haslbeck and Waldorp (2020). Haslbeck, Borsboom, and Waldorp (2018)
introduced Moderated Network Models (MNMs) in more detail for the special case of continuous
variables.

2.3 Applicability of Different Approaches

Table 1 provides an overview of the applicability of the so far discussed methods. The NCT allows
one to compare pairs of GGMs, Ising models and MGMs across two groups. The hypothesis test
based on Fisher’s z-transformation of partial correlations is only applicable to pairs of GGMs. The
FGL allows to compare two or more groups of GGMs. In principle, it could also be implemented
for other models, however, currently no such implementation is available. The method based on

3
the Bayes factor compares the hypotheses that two or more partial correlations are the same or
not. The method based on thresholding the posterior of differences between partial correlations
can compare pairs of GGMs. Using semi-parametric copula models, these two approaches can be
extended to ordinal and binary data, however, they cannot handle nominal categorical variables
with more than two categories. The SEM-based approach using partial pruning can compare two
or more groups and is applicable to the Gaussian and the Ising model. However, detecting group
differences between Ising models is only computationally feasible if the number of variables is
relatively small (up to around 10).

Data Type
Gaussian Ising Mixed >2 Groups
Ś
NCT X X X
Ś Ś Ś
Fisher’s Z X
Ś Ś
FGL X X
˚ Ś
Bayes Factor X X X
X˚
Ś
Post Diff X X
X˚
Ś
Partial Pruning X X
Moderation X X X X

Table 1: An overview of the adaptability of the applicability of the discussed methods; NCT
= Network Comparison Test, Fisher = Hypothesis test based on Fisher’s Z transform; FGL =
Fused Graphical Lasso, Post Diff = Thresholding Posterior of Differences. The asterisk for Partial
Pruning indicates that the method is only feasible for relatively few variables. The asterisk for the
Bayesian methods indicates that these methods cannot handle nominal categorical variables with
more then two categories.

Finally, the moderation approach allows to compare MGMs (and therefore also GGMs and
Ising models) for more than two groups. Note that the NCT, Fisher’s method, and the posterior
difference method can be applied repeatedly to compare more than two groups. However, these
methods are not able to do that in a single shot like the FGL, the Bayes factor method and the
moderation method.

3 Comparing Performance of Different Approaches

In this section we compare the performance of all above discussed approaches in scenarios that are
typical for applied research. The goal is both to compare the relative performance of the methods
and to provide a general overview of how difficult it is to estimate group differences, which is
currently missing in the literature.

3.1 Simulation Setup

I focus on estimating differences in GGMs and Ising models across two groups, which represents
the most common data analysis problems in applied research and allows us to compare all above
discussed approaches. To create group differences in the true models, I use the following approach:
we first create a network model for Group 1, and then change a number of parameters in this
model to obtain the model for Group 2. We then aim to recover those differences.
The initial network of Group 1 is created such that it resembles the network models encountered
in typical psychological applications: I choose network models with p “ 17 variables, because this
is the median number (range 6 ´ 63) of variables in a recent network re-analysis (Haslbeck & Fried,
2017). In the absence of any general knowledge about network structures in psychology we generate
the pairwise dependency structure with a random graph with edge probability P pedgeq “ 0.2. That
is, on average around 27 of the 136 possible edges are present. In the case of the GGM we assign
a draw from a uniform distribution Up´0.2, 0.4q to each of the edges, which results in values
of partial correlations that are commonly observed in empirical data. We set all intercepts to
zero. The interaction parameters of the Ising model are roughly three times larger than partial

4
correlations (see Appendix A), and we therefore take draws from Up´0.6, 1.2q in the case of the
Ising model. To ensure that all binary variables have a large enough variance, we set the threshold
of each node to the negative sum of edge-weights that are connected to it divided by 1.5. This
results in largely negative thresholds, which is what is observed empirically when estimating the
Ising model in the t0, 1u domain.
In the second step, we create the model of Group 2 by randomly picking 20 interaction pa-
rameters of the model of Group 1, and adding ∆θ to each of them. This means that 20 out of
136 possible edges are different across groups, which we consider a reasonable scenario for prac-
tical applications. Note that this is the level of sparsity that really matters in this simulation,
since we are evaluating the performance in estimating group differences. In the GGMs, we vary
∆θ P t0.05, 0.1, 0.2u, in the Ising models we vary ∆θ P t0.15, 0.3, 0.6u. In the Gaussian case, in
rare occasions the partial correlation matrices were not positive definite. In these cases we sampled
repeatedly until the partial correlation matrices of both groups were positive definite.
In order to study performance as a function of sample size, we vary the sample size of each
group n P t20, 37, 68, 233, 431, 795, 1467, 2708, 5000u, which are on a logarithmic scale from 20 to
500. We chose this sequence because it both covers the full range from extremely poor to near-
perfect performance, and represents sample sizes that are typical in psychological applications.
We use the above mentioned methods with a number of different specifications. All methods
include tuning parameters that allow one to make the method more liberal/conservative, and which
exact tuning parameter is selected in a simulation study is therefore somewhat arbitrary. Here we
chose the standard values that are often used in practice to get a rough overview of the performance
of all methods.
We run the NCT as implemented in the R-package NetworkComparisonTest (version 2.2.1; van
Borkulo et al., 2017), where we use 1000 permutations (250 for the Ising model to render the
simulation feasible), use the EBIC for model selection with γ “ 0.25 and evaluate its performance
with significance thresholds α “ 0.05 and α “ 0.01. The FGL is estimated with the R-package
EstimateGroupNetwork (version 0.2.2; Costantini & Epskamp, 2017) using sequential search that
first selects the λ for the graphical lasso penalty, and then the λ for the fused lasso penalty. To select
regularization parameters in the FGL we use 10-fold cross-validation and the EBIC with γ “ 0.25.
We run the Partial Pruning method with α “ 0.05 and α “ 0.01 using the implementation in the
R-package psychonetrics (version 0.7.1; Epskamp, 2020). For the two Bayesian methods we use the
implementations in the BGGM -package (version 2.0.0; Williams & Mulder, 2019). For the Bayes
factor method we define the cut-off value above which a difference is considered reliable as 1. In
this method we set the standard deviation of the Gaussian prior distribution of the group difference
to 0.2. The second Bayesian method uses the posterior of the differences between both groups,
which we threshold to zero if its 95% credible interval overlaps with zero. Finally, we use the
R-package mgm version (version 1.2-9; Haslbeck & Waldorp, 2020) for the implementation of the
moderation approach. We select the regularization parameter with either 10-fold cross-validation
or the EBIC with γ “ 0.25. We also compare the performance of the algorithm with the AND
and the OR-rule. This is especially relevant in the implementation of the mgm-package, since
it runs a nodewise regression on each node, including the moderator variable. This means that
the regression on the moderator variable includes many terms, which renders the AND-rule very
conservative. All simulations were performed with R version 3.6.1.
I evaluate performance with three different measures. First, sensitivity, which is the probability
that a true group difference is recovered. Second, precision, which is the probability that an
estimated group difference is a true group difference. These two measures capture how well the
methods estimate the presence/absence of group differences. However, they provide a poor measure
of how close a given estimate is the true group difference. For example, if the true group difference
is 0 and we estimate it to be 0.0001, this error will impact precision a lot, even through the error
is very small, and in practice it is unlikely that one would interpret such a group difference. We
therefore also consider the estimation error, which we define as the average absolute value of the
difference between true group differences and estimated group differences. For example, if there
are two true group differences ∆θ1,2 “ 0.1, ∆θ6,3 “ 0.3 and the corresponding estimates are equal
to 0 and 0.2, then the average estimation error is equal to p|0.1 ´ 0| ` |0.3 ´ 0.2|q{2 “ 0.1. We
report this estimation error separately for group differences that are present/absent in the true
model.
This design is run for 200 iterations. The code to reproduce all simulation results and Figures
in this paper is available on Github: [Link]

5
3.2 Simulation Results
Figure 1 displays sensitivity, precision, and absolute estimation error for each method for detecting
group differences in GGMs, as a function of the size of the group difference ∆θ and sample size
n. We first focus on the estimation of true group differences. We see that sensitivity increases
for all methods, as expected, however at different rates. The larger the true group difference ∆θ
the quicker sensitivity increases as a function of n, as expected. Within each ∆θ-variation, the
estimation methods differ in how quickly sensitivity increases with n. That is, the methods differ
in how liberal/conservative they are. What stands out is that the MNM with cross-validation (CV)
and the OR-rule is most the liberal, and that the MNM with EBIC and the AND-rule is the most
conservative. Also, the FGL approaches seem to be more liberal than the NCT, which in turn is
more liberal than the Fisher’s method. The latter is followed by the Bayesian methods which have
a similar sensitivity as the partial pruning method. While sensitivity tells us how well a method is
performing in determining whether a group difference is present or absent, it does not tell us how
close the estimate is to the true difference.

NCT, α = 0.05 FGL EBIC, γ = 0.25 MNM EBIC, γ = 0.25 + AND

Fisher, α = 0.05 FGL CV MNM CV + AND
Partial Pruning, α = 0.05 BGGM, Bayes Factor MNM EBIC, γ = 0.25 + OR
BGGM, Post Diff MNM CV + OR

Sensitivity Precision Estimation Error (present) Estimation Error (absent)

1.0 1.0 0.25 0.25

0.8 0.8 0.20 0.20

∆θ = 0.05

0.6 0.6 0.15 0.15

0.4 0.4 0.10 0.10

0.2 0.2 0.05 0.05

0.0 0.0 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000
1.0 1.0 0.25 0.25

0.8 0.8 0.20 0.20

∆θ = 0.10

0.6 0.6 0.15 0.15

0.4 0.4 0.10 0.10

0.2 0.2 0.05 0.05

0.0 0.0 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

1.0 1.0 0.25 0.25

0.8 0.8 0.20 0.20

∆θ = 0.20

0.6 0.6 0.15 0.15

0.4 0.4 0.10 0.10

0.2 0.2 0.05 0.05

0.0 0.0 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

Sample Size per Group Sample Size per Group Sample Size per Group Sample Size per Group

Figure 1: The sensitivity, precision, and estimation errors (separately group differences that are
present or absent in the true model) for estimating differences in interaction parameters across
groups, for the compared approaches, as a function of ∆θ and the number of observations n in
each group. Precision is only displayed for those scenarios in which it was possible to calculate in
at least 90 iterations. This occurs for scenarios with small n and conservative approaches.

The estimation error for true group differences shown in the third column of Figure 1 shows
how close the estimates are to the true group differences. We see that the methods without reg-
ularization have a relatively large estimation error when n is small. The regularized methods, on
the other hand, seem to never have an estimation error that is larger than the true group difference
∆θ. This is because methods without regularization can show errors both due to over- and under-
estimation, while the methods with regularization tend to show only error due to underestimation.

6
Therefore the largest error they can show is ∆θ. While the estimation error provides additional
information to sensitivity, the two measures are of course strongly related: The higher sensitivity,
the lower the estimation error for present group differences.
Next, we consider how well absent group differences are estimated. We only display precision
if edges were estimated to be present in at least 90% of the iterations. This is to avoid precision
estimates that are based only on a few iterations. We also do not show the performance of the
Partial Pruning method for n “ 20, since the underlying model was not always identified in this
scenario. We see that precision is slightly higher for larger group differences, which makes sense
because the presence of true group differences is easier to estimate. We also see that the methods
roughly stack up in the reverse order compared to sensitivity, as one would expect from methods
that offer different trade-offs between sensitivity and precision. The MNM methods with AND-
rule are very conservative, because they also include a regression on the moderator variable, which
includes a very large number of terms and therefore the regularization sets most terms to zero.
While all methods seem to converge to a sensitivity of 1, not all methods seem to converge to a
precision of 1. This is at least the case for the FGL approaches and some of the MNM approaches.
While this is a downside of these methods, precision does not tell us how big those errors are, that
is, how large the estimates of those false positive group differences are. We display the estimation
errors for absent group differences in the fourth column of Figure 1. For low n the unregularized
methods show high estimation errors, which tend to 0 as n increases. The regularized methods
show very low estimation errors for all n. This shows that while these methods incorrectly estimate
some group difference to be present, these incorrect estimates are extremely small. We also see
that this estimation error does not differ across ∆θ which makes sense, because the absent edges
do not vary across these scenarios.
Figure 2 displays sensitivity, precision, and estimation errors for each method for detecting
group differences in the Ising model, as a function of the size of the group difference ∆θ and
sample size n. We again first focus on the estimation of present group differences. Similarly to the
GGM, sensitivity increases fast with n with larger true group differences ∆θ. For the estimation of
group differences in the Ising model only the NCT and the MNM methods are applicable. Similarly
to the GGM-results above, the most liberal method is the MNM with CV and the OR-rule, and
the most conservative method is the MNM with EBIC and the AND-rule, for the same reasons as
in the GGM-case. The remaining methods show comparable sensitivity. The estimation error for
present group differences never exceeds the true group difference, which is because both the NCT
and MNM methods use regularization which biases estimates of group differences towards zero.
Next, we consider the estimation of absent group differences. As expected the precision of the
methods stacks up in reverse order compared to sensitivity. Again similarly to the GGM, some
MNM methods do not seem to approach 1 as n increases, but the estimation errors in the fourth
column show that these errors are very small across all n.
To summarize, we have seen that no method clearly outperforms any other method, but that
they offer different trade-offs between sensitivity and precision. The results above may suggest that
these trade-offs are mostly a function of the type of method. However, there are also considerable
differences within a given method, depending on which tuning parameters are chosen. We see
this for the FGL and MNM methods. In Appendix C we report additional simulation results for
other choices of tuning parameters of the remaining methods (e.g. α “ 0.01 vs α “ 0.05). This
shows that the trade-off between sensitivity/precision is both a function of the type of method and
the particular choice of tuning parameters. Some methods do not seem to converge in precision,
however, the estimation errors show that the errors that lead to low precision are very small. In
general, we saw that a considerable number of observations is necessary in each group in order
to recover small group differences. If group differences are large, however, they can be picked up
even with small sample sizes per group around n “ 100 to 200. Finally, assuming that the relative
scaling of GGM and Ising model parameters I chose is reasonable, it is much harder to detect group
differences in the Ising model.

7
 NCT, α = 0.05 BGGM, Bayes factor MNM EBIC, γ = 0.25 + AND
NCT, α = 0.01 BGGM, Post diff MNM CV + AND
MNM EBIC, γ = 0.25 + OR
MNM CV + OR

Sensitivity Precision Estimation Error (present) Estimation Error (absent)

1.00 1.00 1.00 1.00

0.75 0.75 0.75 0.75

∆θ = 0.15

0.50 0.50 0.50 0.50

0.25 0.25 0.25 0.25

0.00 0.00 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000
1.00 1.00 1.00 1.00

0.75 0.75 0.75 0.75

∆θ = 0.30

0.50 0.50 0.50 0.50

0.25 0.25 0.25 0.25

0.00 0.00 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000
1.00 1.00 1.00 1.00

0.75 0.75 0.75 0.75

∆θ = 0.60

0.50 0.50 0.50 0.50

0.25 0.25 0.25 0.25

0.00 0.00 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000
Sample Size per Group Sample Size per Group Sample Size per Group Sample Size per Group

Figure 2: The sensitivity, precision, and estimation errors (separately group differences that are
present or absent in the true model) for estimating differences in interaction parameters across
groups, for the compared approaches, as a function of ∆θ and the number of observations n in
each group. Precision is only displayed for those scenarios in which it was possible to calculate in
at least 90 iterations. This occurs for scenarios with small n and/or conservative approaches.

4 Tutorial on Estimating Group Differences in Mixed Data

with the Moderation Approach
In this section I show how to use the R-package mgm (Haslbeck & Waldorp, 2020) to estimate
differences across three groups in mixed data using the moderation approach. To be able to share
the data and therefore to make the tutorial fully reproducible for the reader, I use simulated data.
The data consists of four continuous variables, two (nominal) categorical variables with 2 and 3
categories, and a categorical moderator with 3 categories which serves as the grouping variable.
Appendix B describes how I generated the dependencies between these variables.
The data are automatically loaded with the mgm-package and can be accessed in the object
dataGD. The data set contains 3000 observations, 1000 in each of the three groups of the grouping
variable in column 7:

> dim(dataGD)
[1] 3000 7
> head(dataGD)
x1 x2 x3 x4 x5 x6 x7
[1,] 0.214 0.157 0 -0.624 1 0.105 1
[2,] 0.480 0.743 1 -0.273 1 0.428 1
[3,] 0.088 -0.129 1 1.326 0 -0.963 1
[4,] 0.444 -0.487 1 1.670 0 0.641 1

8
[5,] -0.363 -0.080 1 0.306 2 0.376 1
[6,] 0.123 -0.940 0 -0.605 0 1.000 1

In order to detect group differences in the Mixed Graphical Model describing the relationships
between variables X1 , X2 , . . . , X6 we fit a moderated MGM with the grouping variable X7 being
specified as a categorical moderator:

library(mgm)
mgm_obj <- mgm(data = dataGD,
type = c("g", "g", "c", "g", "c", "g", "c"),
level = c(1, 1, 2, 1, 3, 1, 3),
moderators = 7,
lambdaSel = "EBIC",
lambdaGam = 0.25,
ruleReg = "AND")

The argument type indicates the type of variable ("g" for continuous-Gaussian, and "c" for
categorical) and level indicates the number of categories of each variable, which is set to 1 by
default for continuous variables. The moderators argument specifies that the variable in the 7th
column is included as a moderator. Since we specified via the type argument that this variable is
categorical, it will be treated as a categorical moderator. The remaining arguments specify that
the regularization parameters in the `1 -regularized nodewise regression algorithm used by mgm
are selected with the EBIC with a hyperparameter of γ “ 0.25 and that estimates are combined
across nodewise regressions using the AND-rule.
In order to inspect the MGMs in the three groups, we need to condition the moderated MGM
on the values of the three groups. This can be done with the function condition(), which takes
the moderated MGM object and a list specifying on which values of which variables the model
should be conditioned on. Here we only have a single moderator variable (X7 ) and we condition on
each of its values t1, 2, 3u which represent the three groups, and save the three conditional pairwise
MGMs in the list object l_mgm_cond:

l_mgm_cond <- list()

for(g in 1:3) l_mgm_cond[[g]] <- condition(mgm_obj, values = list("7" = g))

We can now inspect the pairwise MGM in each group similar to when fitting a standard pairwise
MGM (for details see Haslbeck & Waldorp, 2020). Here we choose to visualize the strength of
dependencies in the three MGMs in a network using the qgraph package (Epskamp et al., 2012).
We provide the three mgm-objects as an input and set the maximum argument in qgraph() for
each visualization to the maximum parameter across all groups to ensure that the visualizations
are comparable.

library(qgraph)

v_max <- rep(NA, 3)

for(g in 1:3) v_max[g] <- max(l_mgm_cond[[g]]$pairwise$wadj)

par(mfrow=c(1, 3))
for(g in 1:3) {
qgraph(input = l_mgm_cond[[g]]$pairwise$wadj,
[Link] = l_mgm_cond[[g]]$pairwise$edgecolor,
layout = "circle", mar=c(2, 3, 5, 3),
maximum = max(v_max), vsize = 16, esize = 23,
[Link] = TRUE, [Link] = 3)
mtext(text = paste0("Group ", g), line = 2.5)
}

The resulting network visualization is shown in Figure 3. The edges represent conditional
dependence relationships and their width is proportional to their strength. The green (red) edges
indicate positive (negative) linear relationships. The grey edges indicate relationships involving
categorical variables, for which no sign is defined (for details see Haslbeck & Waldorp, 2020). We
see that there are conditional dependencies of equal strength between variables X1 ´ X3 , X3 ´ X4
and X4 ´ X6 in all three groups. However, the linear dependency between X1 ´ X2 differs across

9
groups: it is negative in Group 1, positive in Group 2 and almost absent in Group 3. In addition,
there is no dependency between X3 ´ X5 in Group 1, but there is a dependency in Groups 2
and 3. Note that the equal strength in dependencies between those variables in Groups 2 and 3
does not mean that the exact nature of these dependencies is the same. As with pairwise MGMs,
it is possible to inspect the (non-aggregated) parameter estimates of these interactions with the
function showInteraction().

Group 1 Group 2 Group 3

1 1 1

0.25 0.39 0.06

7 2 7 2 7 2

0.32 0.32 0.32

6 3 6 3 6 3

0.4 0.4 0.4

0.31 0.39 0.31 0.34 0.31

5 4 5 4 5 4

Figure 3: The conditional MGMs in the three groups obtained by conditioning the moderated
MGM on the values of the grouping variables. Green edges indicate positive linear relationships, red
edges indicate negative linear relationships, grey edges indicate relationships involving categorical
variables, for which no sign can be defined. The width of edges is proportional to the strength of
the dependency.

Alternatively to conditioning the moderated MGM one can also inspect the parameters of the
moderated MGM directly. In this case it is important to keep the interpretation of pairwise and
moderation effects in mind. If a given pairwise dependency is moderated (i.e., differs across groups)
then the pairwise interaction indicates the pairwise interaction of the reference group. The smallest
value of the categorical moderator variable is used as the reference category by default.

5 Discussion
In this paper I introduced an alternative way to detect group differences in statistical network
models by including the grouping variable as a moderator in the model. This method is attrac-
tive because it allows to compare networks across more than two groups, and because it is easy
to implement for many models. In addition, I provided a simulation study that evaluated the
performance of the moderation method and existing methods. The results showed that different
methods provide a different trade-off between sensitivity and precision, and that these differences
are determined both by the type of method and the choice of tuning parameters. I also showed that
while some methods perform poorly in precision, their estimation error is actually very low, which
means that incorrectly detected group differences are very small. Given that the scaling between
GGM and Ising parameters chosen here is reasonable, we also saw that it is easier to recover group
differences in the GGM than in the Ising model. Finally, I provided a fully reproducible tutorial
on how to compare an MGM across three groups.
There are several limitations about the reported simulation study that require discussion. First,
the group sizes are equal, which is typically not the case in practice. We assumed that the data is
Gaussian, however in practice data are often skewed. Similarly, in the Ising model we generated
the true models such that all variables have reasonably large variance, which is often not the case in
empirical data. The performance reported in the simulation study should therefore be interpreted
as the best case scenario. Also, we only considered comparisons between pairs of groups, in order to
be able to compare most of the currently used methods. However, we expect that the performance
is similar for three or more groups, as long as the sample size in each group remains constant.
The effect sizes for group differences (or moderation effects) are typically small in observational
behavioral data sets (e.g., Chaplin, 1997; McClelland & Judd, 1993; Sherman & Pashler, 2019).
The true effect sizes are therefore close to ∆θ “ 0.05 or ∆θ “ 0.15 for Gaussian and Ising models,

10
respectively. This suggests that one needs a relatively large sample size to estimate group differences
reliably in observational data. Throughout the simulation study we focused on the comparison of
interaction parameters. However, one could also be interested in mean differences. This is possible
in the NCT, since a sampling distribution is available for any parameter. And it is possible in the
MNM approaches, since the grouping variables is also included as a main effect in the model. To
my knowledge, conditional mean differences can currently not be tested with the FGL approach
and the Bayesian approaches implemented in the BGGM package.
In the implementation of the moderation approach used in this paper, group differences were
estimated regularized GLMs with interaction terms. This means that one does not have a theo-
retically guaranteed false positive rate α or confidence intervals. However, one could also perform
hypothesis tests on the interaction parameters, or one could use the desparsified LASSO (Van de
Geer, Bühlmann, Ritov, Dezeure, et al., 2014) to obtain unbiased sampling distributions to con-
struct confidence intervals.
Interesting avenues for future research would be to extend existing methods to increase their
applicability. For example, one could extend the implementation of the FGL or the Bayesian
methods to also include MGMs. The moderation based method could be improved by using a
more elegant way to integrate estimates across nodewise regressions. Currently, the AND-rule is
extremely conservative, because it includes the regression on the grouping variable, which has a
huge number of parameters. On the other hand, the OR-rule is very liberal, since one nonzero
estimate out of three is sufficient to estimate a group difference to be present. A better aggregation
rule could be to use an AND-rule only on the regressions that do not predict the grouping variable.
To summarize, I introduced a new method to detect group differences between network models
based on moderation analysis, which goes beyond existing methods in that it allows to compare
networks across several groups and grouping variables, and is available for the GGMs, the Ising
models, and MGMs. In addition, I provided the first simulation study comparing the performance
of existing methods to estimate group differences. I hope that these results help applied researchers
to choose the best method in a given situation and to better plan studies that involve investigating
group differences.

Acknowledgements
I would like to thank Claudia van Borkulo, Giulio Costantini, Fabian Dablander, Sacha Epskamp,
Lourens Waldorp, and Donald Willliams for helpful discussions and comments on earlier versions
of this paper.

References
Chaplin, W. F. (1997). Personality, interactive relations, and applied psychology. In Handbook of
personality psychology (pp. 873–890). Elsevier.
Costantini, G., & Epskamp, S. (2017). Estimategroupnetwork: Perform the joint graphical lasso
and selects tuning parameters [Computer software manual]. Retrieved from [Link]
.[Link]/package=EstimateGroupNetwork (R package version 0.1.2)
Danaher, P., Wang, P., & Witten, D. M. (2014). The joint graphical lasso for inverse covari-
ance estimation across multiple classes. Journal of the Royal Statistical Society: Series B
(Statistical Methodology), 76 (2), 373–397.
Epskamp, S. (2020). Psychometric network models from time-series and panel data. Psychometrika,
1–26.
Epskamp, S., Cramer, A. O., Waldorp, L. J., Schmittmann, V. D., Borsboom, D., et al. (2012).
qgraph: Network visualizations of relationships in psychometric data. Journal of Statistical
Software, 48 (4), 1–18.
Epskamp, S., Isvoranu, A.-M., & Cheung, M. (2020). Meta-analytic gaussian network aggregation.
Retrieved from [Link]
Epskamp, S., Maris, G. K., Waldorp, L. J., & Borsboom, D. (2016). Network psychometrics. arXiv
preprint arXiv:1609.02818 .
Epskamp, S., Waldorp, L. J., Mõttus, R., & Borsboom, D. (2018). The gaussian graphical model
in cross-sectional and time-series data. Multivariate Behavioral Research, 53 (4), 453–480.

11
Fisher, R. A. (1915). Frequency distribution of the values of the correlation coefficient in samples
from an indefinitely large population. Biometrika, 10 (4), 507–521.
Foygel, R., & Drton, M. (2010). Extended bayesian information criteria for gaussian graphical
models. In Advances in neural information processing systems (pp. 604–612).
Friedman, J., Hastie, T., & Tibshirani, R. (2008). Sparse inverse covariance estimation with the
graphical lasso. Biostatistics, 9 (3), 432–441.
Fritz, J., Fried, E. I., Goodyer, I. M., Wilkinson, P. O., & van Harmelen, A.-L. (2018). A network
model of resilience factors for adolescents with and without exposure to childhood adversity.
Scientific Reports, 8 (1), 1–13.
Haslbeck, J. M. B., Borsboom, D., & Waldorp, L. (2018). Moderated network models. arXiv
preprint arXiv:1807.02877 .
Haslbeck, J. M. B., & Fried, E. I. (2017). How predictable are symptoms in psychopathological
networks? a reanalysis of 18 published datasets. Psychological Medicine, 47 (16), 2767–2776.
Haslbeck, J. M. B., & Waldorp, L. J. (2020). mgm: Estimating time-varying mixed graphical
models in high-dimensional data. Journal of Statistical Software, 93 (8), 1–46. doi: 10.18637/
jss.v093.i08
McClelland, G. H., & Judd, C. M. (1993). Statistical difficulties of detecting interactions and
moderator effects. Psychological Bulletin, 114 (2), 376.
Meinshausen, N., Bühlmann, P., et al. (2006). High-dimensional graphs and variable selection
with the lasso. The Annals of Statistics, 34 (3), 1436–1462.
Sherman, R., & Pashler, H. (2019). Powerful moderator variables in behavioral science? don’t bet
on them (version 3).
doi: [Link]
van Borkulo, C. D., Borsboom, D., Epskamp, S., Blanken, T. F., Boschloo, L., Schoevers, R. A., &
Waldorp, L. J. (2014). A new method for constructing networks from binary data. Scientific
Reports, 4 (1), 1–10.
van Borkulo, C. D., Boschloo, L., Borsboom, D., Penninx, B. W., Waldorp, L. J., & Schoevers,
R. A. (2015). Association of symptom network structure with the course of depression.
JAMA Psychiatry, 72 (12), 1219–1226.
van Borkulo, C. D., Boschloo, L., Kossakowski, J., Tio, P., Schoevers, R. A., Borsboom, D., &
Waldorp, L. J. (2017). Comparing network structures on three aspects: A permutation test.
Manuscript submitted for publication.
Van de Geer, S., Bühlmann, P., Ritov, Y., Dezeure, R., et al. (2014). On asymptotically optimal
confidence regions and tests for high-dimensional models. The Annals of Statistics, 42 (3),
1166–1202.
van Loo, H., Van Borkulo, C. D., Peterson, R., Fried, E., Aggen, S., Borsboom, D., & Kendler,
K. (2018). Robust symptom networks in recurrent major depression across different levels
of genetic and environmental risk. Journal of Affective Disorders, 227 , 313–322.
Williams, D. R., & Mulder, J. (2019). Bggm: A R package for bayesian gaussian graphical models.
Williams, D. R., Rast, P., Pericchi, L. R., & Mulder, J. (2019). Comparing gaussian graphical
models with the posterior predictive distribution and bayesian model selection.
Yang, E., Baker, Y., Ravikumar, P., Allen, G., & Liu, Z. (2014). Mixed graphical models via
exponential families. In Artificial Intelligence and Statistics (pp. 1042–1050).

12
A Relative Scaling of Parameters in GGM and Ising Model
In order to make the simulation results for the GGM and the Ising model somewhat comparable,
I conducted the following small simulation to determine the relative size of parameters in the two
models. To do so, I simulate from multivariate Gaussian distributions with five dimensions that
have a single nonzero correlation, and I vary the value of this correlation between 0.05 and 0.80.
Then I binarize the data at the median and obtain unbiased estimates for the Ising model in the
t0, 1u domain. The simulation is repeated 50 times. Figure 4 displays the relationship between
correlations in GGMs and parameters in the Ising model.

Relative Scaling of GGM and Ising parameters

3.0

2.5
Parameter in Ising Model

2.0

1.5

1.0

0.5

0.0

0.050 0.175 0.300 0.425 0.550 0.675 0.800

Correlation in GGM

Figure 4: The relative scaling of correlations in GGMs and parameters in the Ising model. The
red line indicates the best fitting regression line.

We see that the relationship can be reasonably well approximated with a linear function. The
parameters of the best fitting regression line are β0 “ ´0.16 and β1 “ 3.3 (red line). We therefore
choose the parameters and group differences in the Ising model to be three times larger than in
the GGM.

13
B Generation of Data used in Tutorial
This appendix describes how we generated the data used in the Tutorial in Section 4. To generate
the data in Group 1 we used the DAG shown in Figure 5 for mixed variables.

X1 „ N p0, 1q

1 X2 „ N p0.5X1 , 1q

1
P pX3 “ 1q “ 1`expt0.5X1 u

2 3 X4 „ N p0.5X3 , 1q

P pX5 “ x5 | X3 “ x3 q
is defined by
4 5
x3 “ 0 x3 “ 1
x5 “ 0 0.166 0.166
˜ ¸
x5 “ 1 0.166 0.166
x5 “ 2 0.166 0.166
6

X6 „ N p0.5X4 , 1q

Figure 5: The definition of the mixed DAG used to generate the data for Group 1.

To generate data for Group 2 I changed the linear effect of X1 on X2 from ´0.5 to 0.5, and I
changed the probability table for the relationship between X3 and X5 to

x3 “ 0 x3 “ 1
x5 “ 0 0.166 0.166
˜ ¸
x5 “ 1 0.100 0.233 .
x5 “ 2 0.233 0.100

In order to generate data for Group 3 I adapted the model of Group 1 by removing the linear
effect of X1 on X2 and by changing the probability table for the relationship between X3 and X5
to

x3 “ 0 x3 “ 1
x5 “ 0 0.266 0.066
˜ ¸
x5 “ 1 0.166 0.166 .
x5 “ 2 0.166 0.166

Thus in the population model there are two dependencies that differ across groups. The de-
pendency between the continuous variables X1 and X2 is negative in Group 1, positive in Group
2 and absent in Group 3. The dependency between the categorical variables X3 and X5 is present
in different forms in Groups 2 and 3, and it is absent in Group 1 (all six cells have the same
probability of 61 « 0.166).
I sampled 1000 observations from each DAG and combined them into the dataset used in
Section 4. The code to generate this dataset with the model described above is available on
Github: [Link]

14
C Additional GGM Simulation Results
Here we present the results of additional variations of the algorithms shown in Figure 1. We
show the performance of the NCT, the Fisher’s method, and the Partial Pruning method also
for α “ 0.01 instead of α “ 0.05. We also show the analytic solution of the Bayesian posterior
difference method.

NCT, α = 0.05 Partial Pruning, α = 0.05 BGGM, Bayes Factor MNM EBIC, γ = 0.25 + AND
NCT, α = 0.01 Partial Pruning, α = 0.01 BGGM, Post Diff (sampling) MNM CV + AND
Fisher, α = 0.05 FGL EBIC, γ = 0.25 BGGM, Post Diff (analytical) MNM EBIC, γ = 0.25 + OR
Fisher, α = 0.01 FGL CV MNM CV + OR

Sensitivity Precision Estimation Error (present) Estimation Error (absent)

1.0 1.0 0.25 0.25

0.8 0.8 0.20 0.20

∆θ = 0.05

0.6 0.6 0.15 0.15

0.4 0.4 0.10 0.10

0.2 0.2 0.05 0.05

0.0 0.0 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000
1.0 1.0 0.25 0.25

0.8 0.8 0.20 0.20

∆θ = 0.10

0.6 0.6 0.15 0.15

0.4 0.4 0.10 0.10

0.2 0.2 0.05 0.05

0.0 0.0 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000
1.0 1.0 0.25 0.25

0.8 0.8 0.20 0.20

∆θ = 0.20

0.6 0.6 0.15 0.15

0.4 0.4 0.10 0.10

0.2 0.2 0.05 0.05

0.0 0.0 0.00 0.00

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

20
37
68
126
233
431
795
1467
2708
5000

Sample Size per Group Sample Size per Group Sample Size per Group Sample Size per Group

Figure 6: The sensitivity, precision, and estimation errors (separately group differences that are
present or absent in the true model) for estimating differences in interaction parameters across
groups, for the compared approaches, as a function of ∆θ and the number of observations n in
each group.

The key observation is that the sensitivity/precision trade-off differs also considerably within
methods, depending on the particular choice of hyperparameters.

15