This document outlines how to use the multilayer.ergm package to fit exponential random graph models (ERGM) in R. The multilayer.ergm package extends functionalities of the ergm package to multilayer network applications by providing a way to count local network configurations, or network motifs, that span more than one network layer. This document assumes you have a working knowledge of ERGMs and multilayer networks as described in Chen (2019).
First, install the package from GitHub. To do so, you will need the devtools package. If you are on a Windows machine, you will also need the standalone Rtools. If you do not already have the network, ergm, and statnet.common packages, installing multilayer.ergm will install them as well.
Install the package by entering the following into your R console. The argument build_opts can be removed entirely if you do not require the package tutorial.
devtools::install_github("tedhchen/multilayer.ergm",
build_opts = c("--no-resave-data", "--no-manual"))Of course, we want to load the package, which will also load the network, ergm, and statnet.common packages.
library(multilayer.ergm)In this exercise, we will be using data from the policy network surrounding legislation over control of the chemical industry in Germany in 1980 (ChemG). This data was used by Leifeld and Schneider (2012) in their analysis of policy communication. To see the data associated with this policy network, we can refer to its help file.
help(chemg_data)From the help file, we can see that there are five network data files. Each of these networks have the same number of nodes, representing the 30 actors in this policy system. We might study these networks independently as monoplex networks, but to understand the entire policy process, which is a system that includes interactions within and across all of these networks, we need to move to the multilayer approach.
The multilayer.ergm package includes utility functions that help create multilayer networks and check network objects to see whether they satisfy the criteria for multilayer networks.
A multilayer network has multiple layers, which is defined by the user. Technically, a network object is compatible with the multilayer.ergm package as long as it has a vertex attribute called layer.mem, which stands for layer membership. This vertex attribute should be a numeric vector that identifies the layer each node belongs to. The following line of code will assign the vertex attribute to the network object.
network%v%'layer.mem' <- layer.idsFor this exercise, let us construct a multiplex network with one layer being a network of scientific communication among the actors (sci) and the other layer being perception of influence among the actors (infrep).
In the matrix representation of a multiplex network, the network layers are on the main diagonal of the block matrix, and the off-diagonal blocks are identity matrices that serve to tie actors across layers. The to.multiplex function is a utility function that assembles networks as layers in a multilayer network, assigns layer ids to the nodes, and sets the off-diagonal matrices as identity matrices.
mnet <- to.multiplex(sci, infrep, output = "network")
mnet## Network attributes:
## vertices = 60
## directed = TRUE
## hyper = FALSE
## loops = FALSE
## multiple = FALSE
## bipartite = FALSE
## total edges= 406
## missing edges= 0
## non-missing edges= 406
##
## Vertex attribute names:
## layer.mem vertex.names
##
## No edge attributesWe can check whether the network object is compatible with multilayer terms manually, or with the check.multilayer function.
mnet%v%'layer.mem'## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2
## [36] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2check.multilayer(mnet)## Layer membership for all nodes are properly specified.Now that we have constructed our multilayer network mnet, we can fit an ERGM. Let us work our way up from a monoplex approach to the scientific communication network.
First, we create the network.
scinet <- network(sci)
scinet## Network attributes:
## vertices = 30
## directed = TRUE
## hyper = FALSE
## loops = FALSE
## multiple = FALSE
## bipartite = FALSE
## total edges= 63
## missing edges= 0
## non-missing edges= 63
##
## Vertex attribute names:
## vertex.names
##
## No edge attributesThen fit an ERGM using the ergm function from the ergm package. Let us include the baseline edge term, a term for how similar a pair of actors’ preferences are (edgecov(prefsim)), and a term for the tendency for actors to reciprocate scientific communication (mutual).
mod.monoplex <- ergm(scinet ~ edges + edgecov(prefsim) + mutual,
control = control.ergm(seed = 206424))(In case anyone is interested 206424 is based on weather statistics in State College, PA when I wrote this vignette: -2°C, 0% precipitation, 64% humidity, and a wind speed of 24km/h.)
Let us look at the output:
summary(mod.monoplex)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: scinet ~ edges + edgecov(prefsim) + mutual
##
## Iterations: 2 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges -2.88849 0.29500 0 -9.792 <1e-04 ***
## edgecov.prefsim -0.01555 0.10156 0 -0.153 0.878
## mutual 2.30823 0.44119 0 5.232 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 1206 on 870 degrees of freedom
## Residual Deviance: 429 on 867 degrees of freedom
##
## AIC: 435 BIC: 449.3 (Smaller is better.)Next, we move to the multiplex with the scientific communication network on one layer, and preceptions of influence on the other layer.
We already created this network earlier, but before we can fit the ERGM, we need to construct a network of sampling constraints that gives the multilayer network its multiplex structure. The to.multiplex function can do this as well. Make sure to specify offzeros = TRUE.
free <- to.multiplex(matrix(1, ncol = 30, nrow = 30), matrix(1, ncol = 30, nrow = 30),
output = "network", offzeros = TRUE)Now we can fit ERGMs of the multiplex network. Let us look at what happens when we fit an ERGM to a multiplex network but do not include any terms that capture dependence across layers.
This within-layer model has the same terms as the previous monoplex model, for both network layers. The multilayer.ergm package includes a set of within-layer terms for multilayer networks. They are usually named after their ergm counterparts, but with _layer appended to the end. They also usually require an argument layer, which indicates which layer the term should be computed on. For a list of network terms that are available for multilayer networks, see ?multilayer_terms.
mod.within <- ergm(mnet
~ edges_layer(layer = 1) + edges_layer(layer = 2)
+ edgecov_layer(prefsim, layer = 1) + edgecov_layer(prefsim, layer = 2)
+ mutual("layer.mem", diff = T),
control = control.ergm(seed = 206424),
constraints = ~fixallbut(free))This is the output:
summary(mod.within)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: mnet ~ edges_layer(layer = 1) + edges_layer(layer = 2) + edgecov_layer(prefsim,
## layer = 1) + edgecov_layer(prefsim, layer = 2) + mutual("layer.mem",
## diff = T)
##
## Iterations: 2 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges_layer.1 -2.88179 0.28836 0 -9.994 < 1e-04 ***
## edges_layer.2 -1.27679 0.18001 0 -7.093 < 1e-04 ***
## edgecov.layer.1.prefsim -0.01759 0.09886 0 -0.178 0.85875
## edgecov.layer.2.prefsim 0.17504 0.06211 0 2.818 0.00483 **
## mutual.same.layer.mem.1 2.30022 0.44410 0 5.179 < 1e-04 ***
## mutual.same.layer.mem.2 0.33416 0.21998 0 1.519 0.12875
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 2412 on 1740 degrees of freedom
## Residual Deviance: 1516 on 1734 degrees of freedom
##
## AIC: 1528 BIC: 1560 (Smaller is better.)Remember that layer 1 is sci and layer 2 is infrep. From the output we see that the coefficients for the scientific layer is effectively the same as those from the monoplex model. (In fact, the only difference is due to randomness in the estimation procedure.) This is to be expected, because while two layers are included, we did not model any dependence between them, meaning that they do not influence each other at all.
Now, let us fit a model that includes cross-layer dependence terms. Specifically, let us include, using duplexdyad, a term for cross-layer reinforcement and a term for cross-layer reciprocity. From the code below, we see that the duplexdyad term takes an argument layers = list(1, 2). It means that the two layers included in the dependence term are layers 1 and 2. This format is common to cross-layer dependence terms in the multilayer.ergm package.
mod.cross <- ergm(mnet
~ edges_layer(layer = 1) + edges_layer(layer = 2)
+ edgecov_layer(prefsim, layer = 1) + edgecov_layer(prefsim, layer = 2)
+ mutual("layer.mem", diff = T)
+ duplexdyad(c("e", "f"), layers = list(1, 2)),
control = control.ergm(seed = 206424),
constraints = ~fixallbut(free))And the output:
summary(mod.cross)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: mnet ~ edges_layer(layer = 1) + edges_layer(layer = 2) + edgecov_layer(prefsim,
## layer = 1) + edgecov_layer(prefsim, layer = 2) + mutual("layer.mem",
## diff = T) + duplexdyad(c("e", "f"), layers = list(1, 2))
##
## Iterations: 2 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges_layer.1 -3.3173 0.3186 0 -10.410 <1e-04 ***
## edges_layer.2 -1.3762 0.1837 0 -7.493 <1e-04 ***
## edgecov.layer.1.prefsim -0.0646 0.1023 0 -0.631 0.5278
## edgecov.layer.2.prefsim 0.1798 0.0633 0 2.841 0.0045 **
## mutual.same.layer.mem.1 2.3703 0.4709 0 5.033 <1e-04 ***
## mutual.same.layer.mem.2 0.3224 0.2210 0 1.459 0.1446
## duplexdyad.e 1.2958 0.2860 0 4.530 <1e-04 ***
## duplexdyad.f -0.1568 0.3133 0 -0.501 0.6167
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 2412 on 1740 degrees of freedom
## Residual Deviance: 1493 on 1732 degrees of freedom
##
## AIC: 1509 BIC: 1553 (Smaller is better.)The duplexdyad.e term in the output is the cross-layer reciprocity effect and the duplexdyad.f term is the cross-layer reinforcement effect. (For more detail on the different configurations in the duplex dyad census, see Fig. 4 of Chen (2019).) As these results show, an actor’s tendency to communicate scientific information to another actor is reinforced by whether they think that other actor is influential and vice versa.