Illustrating the multilayer.ergm Package with a Multiplex Policy Network

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 (2021).

Installation

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)

ChemG Policy Network

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.

Creating a Multilayer Network

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.ids

Policy Multiplex Network

For 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 attributes

We 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 2 2 2
## [39] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

check.multilayer(mnet)

## Layer membership for all nodes are properly specified.

ERGMs for Multilayer Networks

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.

Scientific Communication Monoplex Model

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 attributes

Then 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))

## Warning: 'glpk' selected as the solver, but package 'Rglpk' is not available;
## falling back to 'lpSolveAPI'. This should be fine unless the sample size and/or
## the number of parameters is very big.

(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)

## Call:
## ergm(formula = scinet ~ edges + edgecov(prefsim) + mutual, control = control.ergm(seed = 206424))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                 Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges            -2.8724     0.2802      0 -10.251   <1e-04 ***
## edgecov.prefsim  -0.0190     0.1001      0  -0.190    0.849    
## mutual            2.2969     0.4309      0   5.331   <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. MC Std. Err. = 0.6321)

Next, we move to the multiplex with the scientific communication network on one layer, and preceptions of influence on the other layer.

Policy Multiplex, Within-layer Model

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)

## Call:
## ergm(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), constraints = ~fixallbut(free), 
##     control = control.ergm(seed = 206424))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                         Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges_layer.1           -2.88319    0.29473      0  -9.783  < 1e-04 ***
## edges_layer.2           -1.28077    0.17775      0  -7.206  < 1e-04 ***
## edgecov.layer.1.prefsim -0.01421    0.09982      0  -0.142  0.88683    
## edgecov.layer.2.prefsim  0.17588    0.06013      0   2.925  0.00344 ** 
## mutual.same.layer.mem.1  2.29892    0.44353      0   5.183  < 1e-04 ***
## mutual.same.layer.mem.2  0.33446    0.21596      0   1.549  0.12145    
## ---
## 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: 1561  (Smaller is better. MC Std. Err. = 0.437)

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.

Policy Multiplex, Cross-layer Model

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)

## Call:
## ergm(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)), constraints = ~fixallbut(free), 
##     control = control.ergm(seed = 206424))
## 
## Monte Carlo Maximum Likelihood Results:
## 
##                         Estimate Std. Error MCMC % z value Pr(>|z|)    
## edges_layer.1           -3.28344    0.32925      0  -9.972  < 1e-04 ***
## edges_layer.2           -1.38878    0.18792      0  -7.390  < 1e-04 ***
## edgecov.layer.1.prefsim -0.07683    0.10169      0  -0.756  0.44993    
## edgecov.layer.2.prefsim  0.18571    0.06369      0   2.916  0.00355 ** 
## mutual.same.layer.mem.1  2.36387    0.48700      0   4.854  < 1e-04 ***
## mutual.same.layer.mem.2  0.32128    0.22234      0   1.445  0.14847    
## duplexdyad.e             1.32365    0.29477      0   4.490  < 1e-04 ***
## duplexdyad.f            -0.18291    0.30782      0  -0.594  0.55237    
## ---
## 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. MC Std. Err. = 0.5571)

The duplexdyad.e term in the output is the cross-layer reinforcement effect and the duplexdyad.f term is the cross-layer reciprocity effect. (For more detail on the different configurations in the duplex dyad census, see Fig. 4 of Chen (2021).) 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.