# Bayesian exponential random graph models for populations of networks

Brieuc Lehmann<EMAIL_ADDRESS>[ Simon
<EMAIL_ADDRESS>[  Department of Statistics,
University of Oxford, 29 St Giles’, Oxford, OX1 3LB, United Kingdom  Big Data
Institute, University of Oxford, Old Road Campus, Oxford, OX3 7LF, United
Kingdom,  Department of Psychiatry, University of Cambridge, Cambridge, CB2
OAH, United Kingdom,  MRC Biostatistics Unit, University of Cambridge,
Cambridge, CB2 0SR, United Kingdom

###### Abstract

The collection of data on populations of networks is becoming increasingly
common, where each data point can be seen as a realisation of a network-valued
random variable. A canonical example is that of brain networks: a typical
neuroimaging study collects one or more brain scans across multiple
individuals, each of which can be modelled as a network with nodes
corresponding to distinct brain regions and edges corresponding to structural
or functional connections between these regions. Most statistical network
models, however, were originally proposed to describe a single underlying
relational structure, although recent years have seen a drive to extend these
models to populations of networks. Here, we propose one such extension: a
multilevel framework for populations of networks based on exponential random
graph models. By pooling information across the individual networks, this
framework provides a principled approach to characterise the relational
structure for an entire population. To perform inference, we devise a novel
exchange-within-Gibbs MCMC algorithm that generates samples from the doubly-
intractable posterior. To illustrate our framework, we use it to assess group-
level variations in networks derived from fMRI scans, enabling the inference
of age-related differences in the topological structure of the brain’s
functional connectivity.

exponential random graph models,

Bayesian multilevel modelling,

brain network,

###### keywords:

and

##  1 Introduction

The statistical analysis of network data is becoming increasingly commonplace,
with applications across various disciplines, such as epidemiology, social
science, neuroscience and finance (Kolaczyk, 2009). Over the last four
decades, a number of statistical models for networks have been developed,
including stochastic blockmodels (Holland et al., 1983), latent space models
(Hoff et al., 2002), and the focus of this article, exponential random graph
models (ERGMs; Frank and Strauss (1986)).

An exponential random graph model is a set of parametric statistical
distributions on network data (see Schweinberger et al. (2020) for a recent
review). The aim of the model is to characterise the distribution of a network
in terms of a set of summary statistics. These summary statistics are
typically comprised of topological features of the network, such as the number
of edges and subgraph counts. The summary statistics enter the likelihood via
a weighted sum; the weights are (unknown) model parameters that quantify the
relative influence of the corresponding summary statistic on the overall
network structure and must be inferred from the data. ERGMs are thus a
flexible way in which to describe the global network structure as a function
of network summary statistics.

To date, statistical network models, including ERGMs, have largely focused on
the analysis of a single network. Formally, a network consists of a set of
nodes and a set of edges between these nodes. Let
$\mathcal{N}=\\{1,\dots,N\\}$ be a finite set of nodes, each of which may be
associated with covariates $x_{i}\in\mathcal{X}\subseteq\mathbb{R}^{q}$. An
edge from node $i$ to node $j$ is denoted by $Y_{ij}$, so that the network is
encoded by the adjacency matrix $\boldsymbol{Y}=(Y_{ij})_{i,j\in\mathcal{N}}$.
For our purposes, the set of nodes $\mathcal{N}$ and their covariates
$\boldsymbol{x}=\\{x_{1},\dots,x_{N}\\}$ are considered fixed, while the edges
are considered to be random variables. Denote $\boldsymbol{y}$ to be an
instantiation, or outcome, of the random adjacency matrix $\boldsymbol{Y}$ and
write $\mathbb{P}(\boldsymbol{Y}=\boldsymbol{y}):=\pi(\boldsymbol{y})$ for the
probability that $\boldsymbol{Y}$ takes the value $\boldsymbol{y}$. A
statistical network model specifies a parametrised probability distribution on
the adjacency matrix $\pi(\boldsymbol{y}|\boldsymbol{x},\theta)$ where
$\theta$ is a vector of model parameters.

A population of networks consists of $n>1$ adjacency matrices
$\boldsymbol{Y}^{(1)},\dots,\boldsymbol{Y}^{(n)}$ defined on a common set of
nodes $\mathcal{N}$. We will assume for simplicity that the nodal covariates
are the same across networks, though in principle this not need be the case. A
common example of a population of networks arises in neuroimaging, where a
typical study consists of brain data across a number of participants, each
constituting an individual network. Network analyses of the brain can provide
insight into cognitive function by revealing how distinct brain areas work in
conjunction (Fuster, 2006). These analyses aim to identify salient topological
features of the brain’s connectivity structure that are common across
individuals or that differ across groups.

While one could fit a single model to each individual network separately, it
is not straightforward to combine these individual results into a single
coherent result that is representative of the whole population. An alternative
approach is to construct a group-representative network by, for example,
taking the mean of the edges across the individual networks and applying a
threshold to the resulting weighted network (Achard et al., 2006). These
approaches ignore the individual variability present in the networks and,
moreover, typically do not accurately summarise the topological information
across the individual networks (Ginestet et al., 2011).

A more statistical approach is to treat each individual networks as distinct
statistical units arising from a joint probability distribution
$\pi(\boldsymbol{y}^{(1)},\dots,\boldsymbol{y}^{(n)}|\boldsymbol{x},\theta)$
(Ginestet et al., 2017). Here, we propose a multilevel framework for
populations of networks based on Bayesian exponential random graph models. By
pooling information across the individual networks, this framework provides a
principled approach to characterise the relational structure for an entire
population. Our method can be used to infer group-level differences in the
network structure between sets of networks, which we demonstrate on both
simulated networks and real networks derived from resting-state functional
magnetic resonance imaging (fMRI) scans from an ageing study.

Inference for Bayesian ERGMs is challenging due to the double-intractability
of the ERGM posterior distribution; standard MCMC schemes such as the
Metropolis algorithm are not feasible as it is not possible to evaluate the
acceptance ratio. A common approach is to apply the exchange algorithm (Murray
et al., 2006), which was first employed in the context of Bayesian ERGMs by
Caimo and Friel (2011). To perform inference for our framework for populations
of networks, we propose an exchange-within-Gibbs algorithm that combines the
exchange algorithm with the Gibbs sampler (Geman and Geman, 1984) to produce
samples from the target posterior distribution. The parameterisation of
general multilevel models can play an important role in the overall efficiency
of an MCMC scheme (Gelfand et al., 1995; Papaspiliopoulos et al., 2003, 2007).
To improve the mixing properties of the algorithm, we use an ancillarity-
sufficiency interweaving strategy (ASIS) that interweaves between the centered
and non-centered parameterisations Yu and Meng (2011). To further boost
efficiency, we also employ adaptation of the random-walk proposal parameters
in the algorithm (see e.g. Roberts et al. (1997)).

###  1.1 Related work

### Hierarchical ERGMs

Multilevel networks are networks with a nested hierarchical structure such
that nodes may be grouped into subsets of nodes which may further be grouped
into subset of subsets of nodes, and so on. It is worth emphasising that the
hierarchical nature of a multilevel network corresponds to the grouping of
nodes, as opposed to model parameters as might be typical in a Bayesian
hierarchical model. A population of networks represents a two-level network
such that each subset of nodes correspond to a separate network, with no
connections between distinct subsets (see Figure 1). Wang et al. (2013)
proposed ERGMs for multilevel networks, introducing a range of model
specifications to account for a range of multilevel structures for two-level
networks. Slaughter and Koehly (2016) developed Bayesian hierarchical models
for a single group of networks, incorporating covariate information on each
network. The framework described in this article can be seen as an extension
of this class of models to allow for several groups of networks, allowing for
information to be borrowed both within groups and across groups. Yin et al.
(2020) proposed a mixture of ERGMs to model populations of networks in which
the group membership is unknown. Schweinberger and Handcock (2015) introduced
exponential random graph models with local dependence, providing a general
framework encompassing multilevel networks (and thus populations of networks)
and establishing a central limit theorem for this class of models.

Figure 1: A population of networks (bottom) can be seen as a special case of a
multilevel network (top) in which each subset of nodes contains the same
number of nodes and there are no edges between each subset of nodes.

### ERGMs for brain networks

Exponential random graph models have been applied to resting-state fMRI brain
networks (see Simpson et al. (2011) for an early example). Simpson et al.
(2013) constructed group-representative networks by taking the mean of the
parameter estimates from ERGMs fit to each individual network. Sinke et al.
(2016) constructed group-representative networks directly from individual
diffusion tensor imaging (DTI) brain networks and then fit Bayesian ERGMs to
the resulting group networks. Obando and Fallani (2017) applied ERGMs to
functional connectivity brain networks derived from electroencephalographic
(EEG) signals. In each of these approaches, the networks are fit independently
from each other and, unlike our approach, there is no borrowing of information
across networks.

### Other models for populations of networks

Other statistical network models have recently been extended to handle
populations of networks. Sweet et al. (2013) proposed a general framework of
hierarchical network models (HNMs), which encompasses the model described in
this article. They focus on a hierarchical representation of latent space
models (Hoff et al., 2002) applied to social networks. Sweet et al. (2014)
studied stochastic blockmodel in the HNM framework to infer clusters of nodes
shared across networks. Durante et al. (2017) develop an alternative extension
of the latent space model (Hoff et al., 2002) to populations of networks based
on a low-dimensional mixture model representation. Durante et al. (2018)
applied this model in the context of groups of networks to test for
differences. Mukherjee et al. (2017) used graphons to detect clusters among
multiple networks within a population (as opposed to clusters within
networks). Signorelli and Wit (2020) use a model-based clustering method based
on generalized linear (mixed) models to cluster networks that share certain
network properties of interest.

##  2 Model formulation

###  2.1 Exponential random graph models

The family of exponential random graph models define probability distributions
over the space of networks in terms of sets of summary (or sufficient)
statistics. We will focus on the case of undirected, binary networks, with
$Y_{ij}=Y_{ji}\in\\{0,1\\}$. Let $\mathcal{Y}$ be the range of
$\boldsymbol{Y}$, i.e. the set of all possible outcomes. Let
$s(\boldsymbol{y},\boldsymbol{x})$ denote a vector of $p$ summary statistics,
such that each component is a function
$s_{i}:\mathcal{Y}\times\mathcal{X}^{N}\mapsto\mathbb{R}$.

An ERGM is specified by a particular set of $p$ summary statistics and a map
$\eta:\Theta\mapsto\mathbb{R}^{p}$. The probability mass function of
$\boldsymbol{Y}$ under the corresponding ERGM is given by

$\pi(\boldsymbol{y}|,\boldsymbol{x},\theta)=\dfrac{\exp\left\\{\eta(\theta)^{T}s(\boldsymbol{y},\boldsymbol{x})\right\\}}{Z(\theta)}.$
(1)

Here, $\theta\in\Theta\subseteq\mathbb{R}^{p}$ is a vector of $p$ model
parameters that must be estimated from the data and
$Z(\theta)=\sum_{\boldsymbol{y}^{\prime}\in\mathcal{Y}}\exp\left\\{\eta(\theta)^{T}s(\boldsymbol{y}^{\prime},\boldsymbol{x})\right\\}$
is the normalising constant ensuring the probability mass function sums to
one. Given data, that is, a realisation $\boldsymbol{y}$, the goal is to infer
which values of $\theta$ best correspond to the data under this distribution.
To reduce the notational burden, we will henceforth omit the dependence on the
nodal covariates $\boldsymbol{x}$, considering this to be implicit in the
specification of the probability distribution.

The Bayesian formulation of ERGMs (Koskinen, 2004; Caimo and Friel, 2011)
augments the definition in (1) with a prior distribution $\pi(\theta)$ for the
model parameters. Given an observation $\boldsymbol{y}$ of the network,
inference is performed by analysing the posterior distribution
$\pi(\theta|\boldsymbol{y})$:

$\begin{split}\pi(\theta|\boldsymbol{y})&=\dfrac{\pi(\boldsymbol{y}|\theta)\pi(\theta)}{\pi(\boldsymbol{y})}\\\
&=\dfrac{\exp\left\\{\eta(\theta)^{T}s(\boldsymbol{y},\boldsymbol{x})\right\\}\pi(\theta)}{Z(\theta)\pi(\boldsymbol{y})},\end{split}$
(2)

where
$\pi(\boldsymbol{y})=\int_{\Theta}\pi(\boldsymbol{y}|\theta)\pi(\theta)d\theta$
is the model evidence. The posterior distribution is generally not available
in a closed-form expression. This is due to two properties of the posterior:
first, the (standard) intractability of the model evidence, and second, the
intractability of the likelihood via the intractable normalising constant
$Z(\theta)$. Posterior distributions with these two sources of intractability
are referred to as doubly-intractable. We explore how to sample from this
posterior in Section 3.

###  2.2 A framework for populations of networks

The Bayesian exponential random graph model described above provides a
flexible family of distributions for a single network. Our aim is to extend
this to a model for a population of networks by representing each network as a
separate ERGM within a Bayesian multilevel (or hierarchical) model. By pooling
information across individual networks, this approach allows us to
characterise the distribution of the whole population. Let
$\bm{Y}=(\boldsymbol{Y}^{(1)},\dots,\boldsymbol{Y}^{(n)})$ be a set of $n$
networks. Identify each network $\boldsymbol{Y}^{(i)}$ with its own
individual-level (and vector-valued) ERGM parameter $\theta^{(i)}$. Write
$\bm{\theta}=(\theta^{(1)},\dots,\theta^{(n)})$ for the set of individual-
level parameters.

We model each individual network $\boldsymbol{Y}^{(i)}$ as an exponential
random graph with model parameter $\theta^{(i)}$. Importantly, each individual
ERGM must consist of the same set of $p$ summary statistics $s(\cdot)$. We
assume that, conditional on their respective individual-level parameters, the
$\boldsymbol{Y}^{(i)}$ are independent. We place a population-level prior
distribution $\pi(\bm{\theta}|\phi)$ on the individual-level ERGM parameters.
We again assume conditional independence of the $\theta^{(i)}$ conditional on
the hyperparameter $\phi$, and complete the model with a hyperprior
$\pi(\phi)$ on $\phi$. See Figure 2 for a diagrammatic representation of the
full model.

$\boldsymbol{Y}^{(1)}$$...$$\boldsymbol{Y}^{(n)}$$\theta^{(1)}$$\theta^{(n)}$$...$$\phi$
Figure 2: A diagrammatic representation of the hierarchical framework. Each
network $\boldsymbol{Y}^{(i)}$ is modelled as an exponential random graph with
individual-level parameter $\theta^{(i)}$. In turn, each $\theta_{i}$ is
assumed to come from a common population-level distribution with
hyperparameter $\phi$.

#### Prior specification

The choice of prior distributions for single network Bayesian ERGMs has yet to
be studied in any great detail. The appropriate setting of priors is a
challenging task due to the typically high levels of dependence between
parameters Koskinen et al. (2013). Studies thus far have generally assumed
(flat) multivariate normal prior distributions on the model parameters (Caimo
and Friel, 2011; Sinke et al., 2016; Thiemichen et al., 2016). Through the
choice of the prior distribution $\pi(\phi)$, we can encode additional
information on the population of networks, such as group structure. For
conceptual and computational simplicity, we will also assume multivariate
normal priors, though other prior specifications warrant further
investigation.

Let $g_{i}\in\left\\{1,\dots,J\right\\}$ denote the group membership of
network $i$. On each individual-level parameters, we place a multivariate
normal prior with group-level hyperparameters
$(\mu^{(g_{i})},\boldsymbol{\Sigma}^{(g_{i})}_{\theta})$:

$\theta^{(i)}\sim\mathcal{N}(\mu^{(g_{i})},\boldsymbol{\Sigma}^{(g_{i})}_{\theta}),~{}~{}~{}i=1,\dots,n.$
(3)

We complete the model by specifying (partially) conjugate hyperpriors on the
group-level parameters:

$\begin{split}\mu^{(j)}&\sim\mathcal{N}(\mu^{pop},\boldsymbol{\Sigma}_{\mu})\\\
\boldsymbol{\Sigma}^{(j)}_{\theta}&\sim\mathcal{W}^{-1}(\nu_{\theta},\boldsymbol{\Psi}_{\theta})\\\
\mu^{pop}&\sim\mathcal{N}(\mu_{0},\boldsymbol{\Lambda})\\\
\boldsymbol{\Sigma}_{\mu}&\sim\mathcal{W}^{-1}(\nu_{\mu},\boldsymbol{\Psi}_{\mu})\end{split}$
(4)

This allows information to be pooled across the $J$ groups via the population-
level hyperparameter $\mu^{pop}$ while also allowing for differences in the
group-level means.

##  3 Posterior computation

The double-intractability of the ERGM posterior distribution means that
standard MCMC schemes such as the Metropolis algorithm are not suitable. This
is due to the presence of the intractable normalising constants
$Z(\theta^{(i)}$ in the denominator, rendering calculation of the Metropolis
acceptance rates computationally infeasible. Several methods have been
proposed in recent years to perform Bayesian inference in the presence of
intractable normalising constants (see Park and Haran (2018) for a review). We
focus here on the exchange algorithm (Murray et al., 2006), which was employed
to generate posterior samples for single-network Bayesian ERGMs by Caimo and
Friel (2011). We first recap the exchange algorithm in the context of Bayesian
ERGMs before describing a exchange-within-Gibbs scheme to generate samples
from the joint posterior.

Consider a Metropolis update for a single-network Bayesian ERGM. The
acceptance probability for a proposal $\theta^{\prime}$ from current value
$\theta$ requires evaluation of the ratio $Z(\theta)/Z(\theta^{\prime})$,
which is computationally intractable. The exchange algorithm is an MCMC scheme
designed to circumvent this obstacle. This is achieved by introducing an
auxiliary variable $\boldsymbol{y}^{\prime}\sim\pi(\cdot|\theta^{\prime})$,
i.e. a network drawn from the same exponential random graph model with
parameter $\theta^{\prime}$.

The algorithm targets an augmented posterior

$\pi(\theta,\theta^{\prime},\boldsymbol{y}^{\prime}|\boldsymbol{y})\propto\pi(\theta|\boldsymbol{y})h(\theta^{\prime}|\theta)\pi(\boldsymbol{y}^{\prime}|\theta^{\prime})$
(5)

where $\pi(\theta|\boldsymbol{y})$ is the original (target) posterior,
$h(\theta^{\prime}|\theta)$ is an arbitrary, normalisable proposal function,
and $\pi(\boldsymbol{y}^{\prime}|\theta)$ is the likelihood of the auxiliary
variable. For simplicity, we assume $h(\theta^{\prime}|\theta)$ to be
symmetric. Each of the three terms on the right-hand side of Eq. (5) can be
normalised, so the left-hand side is well-defined as a probability
distribution.

The algorithm proceeds as follows. At each iteration, first perform a Gibbs’
update of $(\theta^{\prime},\boldsymbol{y}^{\prime})$ by drawing
$\theta^{\prime}\sim h(\cdot|\theta)$ followed by
$\boldsymbol{y}^{\prime}\sim\pi(\cdot|\theta^{\prime})$. Next, exchange
$\theta$ and $\theta^{\prime}$ with probability
$\min(1,AR(\theta^{\prime},\theta,\boldsymbol{y},\boldsymbol{y}^{\prime}))$,
where

$\begin{split}AR(\theta^{\prime},\theta,\boldsymbol{y},\boldsymbol{y}^{\prime})&=\dfrac{\pi(\theta^{\prime}|\boldsymbol{y})}{\pi(\theta|\boldsymbol{y})}\cdot\dfrac{\pi(\boldsymbol{y}^{\prime}|\theta)}{\pi(\boldsymbol{y}^{\prime}|\theta^{\prime})}\\\
&=\dfrac{\exp\left\\{\theta^{\prime
T}s(\boldsymbol{y})\right\\}\pi(\theta^{\prime})}{\exp\left\\{\theta^{T}s(\boldsymbol{y})\right\\}\pi(\theta)}\dfrac{Z(\theta)}{Z(\theta^{\prime})}\cdot\dfrac{\exp\left\\{\theta^{T}s(\boldsymbol{y}^{\prime})\right\\}}{\exp\left\\{\theta^{\prime
T}s(\boldsymbol{y}^{\prime})\right\\}}\dfrac{Z(\theta^{\prime})}{Z(\theta)}\\\
&=\exp\left\\{[\theta^{\prime}-\theta]^{T}[s(\boldsymbol{y})-s(\boldsymbol{y}^{\prime})]\right\\}\frac{\pi(\theta^{\prime})}{\pi(\theta)}\end{split}$
(6)

Crucially, the ratio of intractable normalising constants cancel out, and so
this acceptance ratio can indeed be evaluated. The stationary distribution of
the Markov chain constructed through this scheme is
$\pi(\theta,\theta^{\prime},\boldsymbol{y}^{\prime}|\boldsymbol{y})$ (Murray
et al., 2006). Thus, by marginalising out $\theta^{\prime}$ and
$\boldsymbol{y}^{\prime}$, the algorithm yields samples from the desired
posterior, namely $\pi(\theta|\boldsymbol{y})$.

Algorithm 1 The exchange algorithm for a Bayesian ERGM (Caimo and Friel, 2011)

number of MCMC iterations $K$, initial value $\theta_{0}$

for $k=1,\dots,K$ do

\- draw $\theta^{\prime}\sim h(\cdot|\theta_{k-1})$

\- draw $\boldsymbol{y}^{\prime}\sim\pi(\cdot|\theta^{\prime})$

\- set $\theta_{k}=\theta^{\prime}$ with probability
$\min\left(1,AR(\theta^{\prime},\theta_{k-1},\boldsymbol{y},\boldsymbol{y}^{\prime})\right)$
$\triangleright$ See Eq. (6)

\- else, set $\theta_{k}=\theta_{k-1}$

end for

At each iteration, the exchange algorithm requires a sample
$\boldsymbol{y}^{\prime}$ from the ERGM $\pi(\cdot|\theta^{\prime})$ in order
to compute the acceptance ratio. Although perfect sampling for ERGMs is
possible, it is computationally impractical except for a few special cases
(Butts, 2018). A pragmatic alternative, employed in Caimo and Friel (2011) and
Wang and Atchadé (2014), is to use the final iteration of a Metropolis-
Hastings algorithm as an approximate sample from $\pi(\cdot|\theta^{\prime})$
(Hastings, 1970; Hunter et al., 2008a). A theoretical justification of this
approach is given by Everitt (2012): under certain conditions, despite using
an approximate sample, the algorithm nevertheless targets an approximation to
the correct posterior distribution. Further, this approximation improves as
the number of iterations of the inner MCMC increases.

###  3.1 The exchange-within-Gibbs algorithm

We now extend the exchange algorithm in order to generate samples from our
full posterior on a population of networks. As the name suggests, the
exchange-within-Gibbs algorithm combines the exchange algorithm with the Gibbs
sampler (Geman and Geman, 1984) to produce samples from the desired posterior.
Note that we can treat the unknown parameters of the model
$(\boldsymbol{\theta},\phi)$ as components of a single multi-dimensional
parameter. We iteratively sample each component from its conditional
distribution given the remaining components.

The full exchange-within-Gibbs scheme is outlined in Algorithm 2. Since each
step samples from the respective full conditional distribution, the algorithm
ensures that the stationary distribution of the resulting Markov chain is
indeed the joint posterior $\pi(\boldsymbol{\theta},\phi|\bm{y})$ (Tierney,
1994). As with the exchange algorithm for the single-network Bayesian ERGM,
the most computationally expensive step is sampling $\boldsymbol{y}^{\prime}$
from $\pi(\cdot|\theta^{\prime})$, i.e. simulating an exponential random graph
with parameter $\theta^{\prime}$. Moreover, this step must be performed for
each of the individual-level parameter $\theta^{(i)}$ updates. Thus, the
computational cost of each iteration increases linearly with the number of
networks in the data. However, since (conditional on $\phi$) the
$\theta^{(i)}$ are independent, these updates may be performed in parallel.
Therefore, with access to a sufficient number of computing cores, the actual
computational time per iteration typically increases sub-linearly with the
number of networks.

Algorithm 2 The exchange-within-Gibbs algorithm for a multilevel Bayesian ERGM

number of MCMC iterations $K$, initial values
$(\phi_{0},\theta^{(1)}_{0},\dots,\theta^{(n)}_{0})$

for $k=1,\dots,K$ do

\- draw $\phi_{k}\sim\pi(\cdot|\boldsymbol{\theta}_{k-1})$

for $i=1,\dots,n$ do

\- update $\theta^{(i)}_{k}$ via the exchange algorithm

end for

end for

#### Choice of parametrisation: centering vs. non-centering

The parametrisation of general multilevel models in the context of MCMC
computation has been studied in some detail (Gelfand et al., 1995;
Papaspiliopoulos et al., 2003, 2007; Yu and Meng, 2011). Here, we discuss the
two most commonly used parametrisations: the ‘centered’ and the ‘non-
centered’. The parametrisation presented thus far is known as the centered
parametrisation (CP) Gelfand et al. (1995); Papaspiliopoulos et al. (2007), in
which the group-level parameters ($\mu,\boldsymbol{\Sigma}_{\theta}$) are
independent of the data $\bm{Y}$:

$\begin{split}\boldsymbol{Y}^{(i)}&\sim
ERGM(\theta^{(i)}),~{}~{}~{}i=1,\dots,n\\\
\theta^{(i)}&\sim\mathcal{N}(\mu,\boldsymbol{\Sigma}_{\theta}),~{}~{}~{}i=1,\dots,n.\end{split}$
(7)

In contrast, the non-centred parametrisation (NCP) can be written as follows:

$\begin{split}\boldsymbol{Y}^{(i)}&\sim
ERGM(\mu+\tilde{\theta}^{(i)}),~{}~{}~{}i=1,\dots,n\\\
\tilde{\theta}^{(i)}&\sim\mathcal{N}(0,\boldsymbol{\Sigma}_{\theta}),~{}~{}~{}i=1,\dots,n.\end{split}$
(8)

The transformation $\theta^{(i)}=\mu+\tilde{\theta}^{(i)}$ confirms the
equivalence of the two parametrisations. Note that the group-level parameter
$\mu$ enters the likelihood directly in (8), so the conditional distribution
of $\mu$ given the remaining parameters has an intractable normalising
constant. As above, this can be dealt with via an exchange update, in this
case requiring simulation of $n$ networks for the normalising constants to
cancel in the acceptance ratio.

The centred parametrisation and the non-centred parametrisation tend to be
complementary: when one performs poorly, the other tends to perform better
(Papaspiliopoulos et al., 2007). However, when the parameters of interest are
the group-level parameters $(\mu,\boldsymbol{\Sigma}_{\theta})$, it is
possible to combine both approaches using an ancillarity-sufficiency
interweaving strategy (ASIS; Yu and Meng (2011)). ASIS works by combining the
updating schemes of the CP and NCP approaches. The ASIS algorithm for a
multilevel Bayesian ERGM is described in Algorithm 3.

Algorithm 3 The ancillarity-sufficiency interweaving strategy (ASIS) algorithm
for a multilevel Bayesian ERGM

number of MCMC iterations $N$, initial values
$(\mu_{0},\tilde{\theta}^{(1)}_{0},\dots,\tilde{\theta}^{(n)}_{0})$

for $k=1,\dots,K$ do

Draw
$\boldsymbol{\Sigma}_{\theta,k}\sim\pi(\cdot|\mu_{k-1},\tilde{\boldsymbol{\theta}}_{k-1})$

for $i=1,\dots,n$ do

Update $\theta^{(i)}_{k-0.5}$ using the exchange algorithm

end for

Draw centered
$\mu_{k-0.5}\sim\pi(\cdot|\boldsymbol{\theta}_{k-0.5},\boldsymbol{\Sigma}_{\theta,k})$

for $i=1,\dots,n$ do

Set $\tilde{\theta}^{(i)}_{k}=\theta_{k-0.5}^{(i)}-\mu_{k-0.5}$

end for

Update non-centered $\mu_{k}$ (see (8)) via the exchange algorithm

end for

#### Proposal adaptation

We use multivariate normal random walk proposals in the respective exchange
updates of both $\theta^{(i)}$ and $\mu$, for example

$h(\theta^{\prime}|\theta_{k-1})=\mathcal{N}(\theta_{k-1},\Sigma)$ (9)

The choice of the the proposal covariance matrix $\Sigma$ is crucial to the
overall efficiency of the MCMC algorithm; we wish to make large proposals that
are likely to be accepted in order to explore the posterior in as few
iterations as possible. A common approach to tuning covariance proposals for a
wide range of random walk based algorithms, including Metropolis-within-Gibbs,
is to target an acceptance rate close to 0.234, with acceptance rates between
0.1 and 0.5 often yielding satisfactory results (Roberts et al., 1997, 2001;
Roberts and Rosenthal, 2009). Since manual tuning of the $n+1$ proposal
covariance matrices would be impractical, we instead implement an adaptive
proposal scheme.

For each proposal, we use a version of the adaptive Metropolis algorithm
(Haario et al., 2001) considered by Roberts and Rosenthal (2009).
Specifically, for the first 1000 iterations, we adapt every 20 iterations,
with proposals of the form

$h_{k}(\theta^{\prime}|\theta_{k-1})=(1-\beta)N\left(\theta_{k-1},(2.38)^{2}\delta_{k}\Sigma_{k}/p\right)+\beta
N\left(\theta_{k-1},(0.1)^{2}\delta_{k}I_{p}/p\right),$ (10)

where $\Sigma_{k}$ is the sample covariance matrix of the posterior samples
$(\theta_{1},\dots,\theta_{k-1})$ and $\delta_{k}$ is an additional scaling
factor that is varied to control the magnitude of the proposals. Following
Roberts and Rosenthal (2009), we set $\beta=0.05$. The role of $\Sigma_{k}$ is
to adapt the direction of the proposals to the MCMC run so far, while
$\delta_{k}$ serves to target an acceptance rate of 0.234. Specifically, we
start with $\delta_{1}=1$ and increase (resp. decrease) $\log(\delta_{k})$ by
$\min(0.5,1/\sqrt{(}k)$ if the acceptance rate was below (resp. above) 0.234
in the previous 20 iterations.

###  3.2 Posterior predictive assessment

Having produced a sufficient number of samples from the posterior
distribution, we then assess whether the model adequately describes the data.
Since determining the distribution of appropriate test quantities is
difficult, assessing such goodness-of-fit for ERGMs is typically performed
graphically (Hunter et al., 2008b). For a single ERGM fit, one can simulate a
large number of networks from the fitted model and compare these ‘posterior
predictive networks’ to the observed network. This comparison is usually done
via a set of network metrics. If a model fits the data well then the network
metrics of the posterior predictive networks should be similar to those of the
observed network.

For a population of networks, we can apply the same principles. To do so, we
choose uniformly at random $S$ values from the posterior samples of the group-
level mean parameters. For each value, we simulate a network from
$\pi(\cdot|\mu^{(s)})$. We can then compare these posterior predictive
networks to the observed networks based on a set of network metrics. For this
purpose, we will use three important network metric distributions that are not
explicitly modelled, namely degree distribution, geodesic distance
distribution (length of shortest paths) and edge-wise shared partners
distribution.

##  4 Results

To illustrate our method, we apply it to a set of simulated networks,
demonstrating that it is capable of recovering the ground truth. We also apply
our method to resting-state fMRI networks from the Cam-CAN project, a study on
healthy ageing (Shafto et al., 2014), finding that differences in network
structure between a group of young individuals and and a group of old
individuals can largely be explained by a reduction in inter-hemispheric
connections in the old. The R scripts used to generate these results can be
found at https://github.com/brieuclehmann/multibergm-scripts.

###  4.1 Simulation

We generated sets of 30-node networks with nodes split into two ‘hemispheres’
of 15 nodes each. We simulated the networks from an exponential random graph
model with three terms: total number of edges (‘edges’), total number of edges
between nodes in the same hemisphere (‘nodematch.hemisphere’), and the
geometrically-weighted edgewise-shared partner (GWESP) statistic
(‘gwesp.fixed.0.9’). The GWESP statistic of a network $y$ is a measure of
clustering and is given by:

$GWESP(\boldsymbol{y})=e^{\tau}\sum_{w=1}^{N}\\{1-(1-e^{-\tau})^{w}\\}EP_{w}(\boldsymbol{y}),$
(11)

where $EP_{w}(y)$ is the number of connected node pairs having exactly $w$
shared partners and $\tau$ is a decay parameter, which we fix at $\tau=0.9$.
The decay parameter attenuates the effect of the number of higher-order
edgewise shared partners relative to lower-order edgewise shared partners.

To illustrate the method, we first simulate networks for a single group with
$n=10,20,50$ networks, before moving onto an example with two groups, such
that $n=20,40,100$ and $n_{1}=n_{2}$, i.e. both groups have the same number of
networks.

#### Single group

To simulate the networks, we first generated individual-level parameters
$\theta_{i}\sim\mathcal{N}(\mu,\Sigma),~{}i=1,\dots,n$ where

$\displaystyle\mu$ $\displaystyle=(-3,0.5,0.5)^{T}$ (12) $\displaystyle\Sigma$
$\displaystyle=\frac{1}{50}\begin{pmatrix}1&-0.5&0\\\ -0.5&0.5&0\\\
0&0&0.5\end{pmatrix}.$ (13)

We then used the ergm R package (Hunter et al., 2008a) to simulate $n$
networks $\boldsymbol{y}_{i}\sim p(\cdot|\theta_{i}),~{}i=1,\dots,n$. The
simulation procedure is based on an MCMC algorithm, initialised at a network
with the prescribed number of nodes and covariates (in this case, hemisphere
labels). With these simulated networks, we applied our exchange-within-Gibbs
algorithm with ASIS (Algorithm 3) to generate 12,000 posterior samples,
adapting the random-walk proposals for the first 1,000 iterations, and
discarding the first 2,000 as burn-in.

Figure 3 displays summaries of the posterior samples for the group-level mean
parameter $\mu$ of the model fit to $n=10$ networks. The true value of $\mu$
is covered by the posterior density, while the trace and autocorrelation plots
indicate that the MCMC has mixed well. To assess the goodness-of-fit, we
generated $S=100$ networks from the model at posterior samples of $\mu$ chosen
uniformly at random. Figure 4 shows the degree distribution, geodesic distance
distribution and edgewise shared partner distribution of these simulated
networks against those to which the model was fit.

Figure 3: Posterior samples produced by the exchange-within-Gibbs algorithm
for the group-level mean parameter, $\mu$, of a single group of ten simulated
networks. The true value of $\mu$ is indicated by the red line. Figure 4:
Graphical goodness-of-fit assessment for a single group of ten simulated
networks. The box plots correspond to the simulated networks, while the
ribbons represent 90% credible intervals corresponding to the posterior
predictive networks. Note that a geodesic distance of infinity between two
nodes means that there is no path connecting the nodes.

To complete our analysis of a single group of networks, we compare the density
of the posterior samples between groups of size $n=10,20,50$. Figure 5
illustrates how the posterior samples of $\mu$ concentrates around the true
value as the number of networks in the group increases.

Figure 5: Posterior density plots for the group-level mean parameters with
varying number of networks. As the number of networks increases, the posterior
concentrates around the true value, depicted by the red vertical line.

#### Two groups

We now consider a multilevel setting in which the networks are split into two
distinct groups $j=1,2$, each with their own group-level mean parameters
$\mu^{(j)}$ but with common group-level covariance $\Sigma$. As above, we
first generated individual-level parameters
$\theta_{i}\sim\mathcal{N}(\mu^{(g_{i})},\Sigma),~{}i=1,\dots,n$, where
$g_{i}\in\\{1,2\\}$ denotes the group membership of the $i^{th}$ network, and
then simulated networks $\boldsymbol{y}_{i}\sim
p(\cdot|\theta_{i}),~{}i=1,\dots,n$. We considered a range of numbers of
networks, $n=10,20,50$ each of the two groups. The true values were

$\displaystyle\mu^{(1)}$ $\displaystyle=(-3,0.5,0.5)^{T}$ (14)
$\displaystyle\mu^{(2)}$ $\displaystyle=(-2.6,0.5,0.2)^{T}$ (15)
$\displaystyle\Sigma$ $\displaystyle=\frac{1}{50}\begin{pmatrix}1&-0.5&0\\\
-0.5&0.5&0\\\ 0&0&0.5\end{pmatrix}.$ (16)

Figure 6 shows the density of the posterior samples for the group-level
parameters $(\mu^{(1)},\mu^{(2)})$ for increasing number of networks $n$ per
group. We see that, as in the single-group setting, the posteriors concentrate
around the true values for each group as the number of networks increases.

Figure 6: Posterior density plots for a two-group model with $n=10,20,50$
simulated networks in each group.

###  4.2 Application to human functional connectivity brain networks

We now turn our attention to a real data example: networks derived from
resting-state fMRI scans of human brains from the Cambridge Centre for Ageing
and Neuroscience (Cam-CAN) research project Shafto et al. (2014), a study on
the effect of healthy ageing on cognitive and brain function. The Cam-CAN
dataset consists of a range of cognitive tests and functional neuroimaging
experiments for approximately 650 healthy individuals aged 18-87. Our aim will
be to compare the functional connectivity structure between the 100 youngest
individuals, aged 18–33, and the 100 oldest individuals, aged 74–87.

Full details of data collection and preprocessing can be found in Lehmann et
al. (2021). To summarise, both structural (T1 and T2) and eyes-closed,
resting-state fMRI scans (261 volumes, lasting 8min 40s) were acquired for
each individual. The fMRI scans were motion-corrected and co-registered to the
respective structural scans and then mapped to the common Montreal
Neurological Institute (MNI) template to ensure comparability across
individuals. The fMRI time series were then extracted from 90 cortical and
subcortical regions of interest (ROIs) from the AAL atlas (Tzourio-Mazoyer et
al., 2002) and adjusted for various confounds using the optimised pipeline of
Geerligs et al. (2017).

To construct networks for each individual, we followed a thresholded
correlation matrix approach. For individual $i$, we computed the pairwise
Pearson correlation between each of the $N=90$ preprocessed time series,
yielding a $N\times N$ correlation matrix $\boldsymbol{C}^{(i)}$. We then
applied a threshold $r$ to $\boldsymbol{C}^{(i)}$ to produce an $N\times N$
adjacency matrix, $\boldsymbol{A}^{(i)}$, with entries:

$\boldsymbol{A}^{(i)}_{kl}=\begin{cases}1&\text{if
}\boldsymbol{C}^{(i)}_{kl}\geq r\quad k,l=1,\dots,N\\\
0&\text{otherwise.}\end{cases}$ (17)

The adjacency matrix defines an individual’s network, $\boldsymbol{y}^{(i)}$,
with an edge between nodes $k$ and $l$ if and only if
$\boldsymbol{A}^{(i)}_{kl}=1$. The threshold $r$ was chosen to yield an
average node degree of 3 across all the networks, as recommended by Fallani et
al. (2017).

We model the resulting population of $n=200$ networks using the framework
described in Section 2.2 with an exponential random graph model with four
terms: total number of edges (‘edges’), total number of edges between nodes in
the same hemisphere (‘nodematch.hemisphere’), total number of edges between
homotopic nodes (mirror ROIs in each hemisphere; ‘nodematch.homotopy’) and the
geometrically-weighted edgewise-shared partner (GWESP) statistic with decay
parameter $tau=0.9$ (‘gwesp.fixed.0.9’). We again assume a common group-level
covariance structure and focus on the group-level mean parameters. The full
model is given by:

$\begin{split}\boldsymbol{y}^{(i)}&\sim
p(\cdot|\theta^{(i)}),~{}i=1,\dots,200,\\\
\theta^{(i)}&\sim\mathcal{N}(\mu^{(g_{i})},\Sigma_{\theta}),~{}i=1,\dots,200,\\\
\mu^{(j)}&\sim\mathcal{N}(\mu^{pop},\boldsymbol{\Sigma}_{\mu}),~{}j=1,2\\\
\boldsymbol{\Sigma}_{\theta}&\sim\mathcal{W}^{-1}(\nu_{\theta},\boldsymbol{\Psi}_{\theta})\\\
\mu^{pop}&\sim\mathcal{N}(\mu_{0},\boldsymbol{\Lambda})\\\
\boldsymbol{\Sigma}_{\mu}&\sim\mathcal{W}^{-1}(\nu_{\mu},\boldsymbol{\Psi}_{\mu}).\end{split}$
(18)

We used the exchange-within-Gibbs algorithm with ASIS to generate 22,000
posterior samples, discarding the first 2,000 samples as burn-in. Figure 7
shows summaries of the posterior samples for $(\mu^{(1)},\mu^{(2)})$, with the
trace and autocorrelation plots demonstrating that the MCMC has mixed well.
The posterior density plots show that the clearest difference between the old
group and the young group was the difference in the parameter associated with
the number of edges between homotomic nodes (‘nodematch.homotopy’). While this
parameter is large and positive for both groups, it is moderately smaller in
the old group, indicating that the propensity for homotopic connections is
lower in old age. On the other hand, there is no clear evidence for group
differences in the remaining parameters. The edges parameters are large and
negative, pointing to the overal sparsity of the networks; the intrahemisphere
parameters (‘nodematch.hemisphere’) are small and positive, indicating a
moderate propensity for connections between nodes in the same half of the
brain; and the GWESP parameters are also positive, indicating a propensity to
form triangles and thus a degree of functional segregation (Bullmore and
Sporns, 2009).

Figure 7: MCMC output from the exchange-within-Gibbs algorithm for the group-
level mean parameters of a population of resting-state fMRI networks from a
group of 100 young individuals and a group of 100 old individuals.

To assess goodness-of-fit, for both groups we generated $S=100$ networks from
the model at posterior samples of $\mu^{(j)}$ chosen uniformly at random.
Figure 8 indicates a reasonable fit for both groups, with the geodesic
distance and edgewise shared partner distributions showing a good
correspondence between the simulated networks and the observed networks. There
appears to be a slight discrepancy in the degree distributions, with the
simulated networks in the young group in particular having fewer nodes of
degree 4 to 6 relative to the observed networks.

Figure 8: Graphical goodness-of-fit assessment for resting-state fMRI networks
from a young group and an old group, fitted in a joint model. The box plots
correspond to the observed networks, while the ribbons represent 95% credible
intervals corresponding to the posterior predictive networks. Note that a
geodesic distance of infinity between two nodes means that there is no path
connecting the nodes.

##  5 Discussion

The main contribution of this article is to introduce a multilevel framework
for modelling populations of networks, along with a novel MCMC procedure for
performing inference with the framework. While the framework itself is a
natural multilevel extension of a single ERGM, the inference procedure is more
involved due to the intractability of the ERGM likelihood and the challenges
associated with MCMC for hierarchical models. We have presented how our
framework can be applied to functional connectivity networks from two groups
of individuals. The framework allows the pooling of information across
individuals within the same group and, depending on the parametrisation,
across groups. By analysing the posterior samples produced by the MCMC
algorithm, we can determine whether there are statistically significant
differences in the connectivity structure between the two groups. Although we
chose here to focus on networks constructed from resting-state fMRI scans, our
framework could also be applied to networks derived from other neuroimaging
modalities such as magnetoencephalography (MEG) or diffusion tensor imaging
(DTI).

The flexibility of the multilevel framework could be exploited to model more
complex group structures such as multiple groups or factorial designs.
Similarly, given multiple network measurements per individual, it is
straightforward (conceptually, if not computationally) to extend the framework
by adding another layer to the model. One could also incorporate individual-
level covariate information into the model. In the context of the healthy
ageing example considered above, this would avoid splitting individuals into
(somewhat arbitrary) groups and provide a mechanism to measure how
connectivity structure varies throughout the lifespan.

Another important extension would be to use weighted exponential random graph
models (Krivitsky, 2012; Desmarais and Cranmer, 2012). These are an extension
of the binary ERGM that can be applied to weighted networks, thus avoiding the
thresholding step in the construction of functional connectivity networks.
Indeed, one version of a weighted ERGM, the generalised exponential random
graph model (GERGM) (Desmarais and Cranmer, 2012) was recently applied to a
20-node functional connectivity network (Stillman et al., 2017). This approach
has the additional advantage of modelling the mean connectivity directly and
thus would avoid any confounding due to differences in mean connectivity.
However, the GERGM is at present extremely computational intensive, rendering
it infeasible for a population of networks.

One of the key challenges in applying our framework to real data is the choice
of which network summary statistics to include in the model. A fully Bayesian
model selection method based on reversible-jump MCMC has been developed for
exponential random graph models on single networks (Caimo and Friel, 2013). A
similar approach could be developed for our framework, though the
computational cost is likely to be prohibitive. A more pragmatic approach
would be to develop a graphical goodness-of-fit method by comparing the
posterior predictive distributions under different models.

The computational cost of our MCMC algorithm is considerable. Using a 20-core
computing cluster (Intel Skylake, 2.6 GHz), the algorithm took 29 hours to
produce the 22,000 posterior samples in the real data example presented above.
The main computational bottleneck lies in simulating the exponential random
graphs at each MCMC iteration. While the computational cost should increase
roughly linearly in the number of networks, Krivitsky and Handcock (2014)
provide empirical evidence indicating that the cost may grow on the order of
$p(N+E)\log(E)$ where $p$ is the number of summary statistics, $N$ is the
number of nodes, and $E$ is the number of edges. It may be possible to reduce
the number of ERGM simulations at each MCMC iteration using noisy Monte Carlo
methods (Alquier et al., 2016). Another promising avenue is variational
inference for ERGMs (Tan and Friel, 2020), which could be extended to our
framework to yield approximate Bayesian inference at a much reduced
computational cost relative to MCMC.

## Acknowledgements

B.L. and S.W. were supported by the UK Medical Research Council [Programme
number U105292687]. B.L. was also supported by the UK Engineering and Physical
Sciences Research Council through the Bayes4Health programme [Grant number
EP/R018561/1] and gratefully acknowledges funding from Jesus College, Oxford.
This research was supported by the NIHR Cambridge Biomedical Research Centre
(BRC-1215-20014). The computational aspects of this research were supported by
the Wellcome Trust Core Award Grant Number 203141/Z/16/Z and the NIHR Oxford
BRC. The views expressed are those of the authors and not necessarily those of
the NHS, the NIHR or the Department of Health and Social Care.

## References

  * Achard et al. (2006) Achard, S., Salvador, R., Whitcher, B., Suckling, J., and Bullmore, E. (2006).  “A Resilient, Low-Frequency, Small-World Human Brain Functional Network with Highly Connected Association Cortical Hubs.”  Journal of Neuroscience, 26(1): 63–72.
URL http://www.jneurosci.org/content/26/1/63

  * Alquier et al. (2016) Alquier, P., Friel, N., Everitt, R., and Boland, A. (2016).  “Noisy Monte Carlo: convergence of Markov chains with approximate transition kernels.”  Statistics and Computing, 26(1): 29–47.
URL https://doi.org/10.1007/s11222-014-9521-x

  * Bullmore and Sporns (2009) Bullmore, E. and Sporns, O. (2009).  “Complex brain networks: graph theoretical analysis of structural and functional systems.”  Nat Rev Neurosci, 10(3): 186–198.
URL http://dx.doi.org/10.1038/nrn2575

  * Butts (2018) Butts, C. T. (2018).  “A perfect sampling method for exponential family random graph models.”  The Journal of Mathematical Sociology, 42(1): 17–36.
  * Caimo and Friel (2011) Caimo, A. and Friel, N. (2011).  “Bayesian inference for exponential random graph models.”  Social Networks, 33(1): 41 – 55.
URL http://www.sciencedirect.com/science/article/pii/S0378873310000493

  * Caimo and Friel (2013) — (2013).  “Bayesian model selection for exponential random graph models.”  Social Networks, 35(1): 11 – 24.
URL http://www.sciencedirect.com/science/article/pii/S0378873312000573

  * Desmarais and Cranmer (2012) Desmarais, B. A. and Cranmer, S. J. (2012).  “Statistical inference for valued-edge networks: The generalized exponential random graph model.”  PloS one, 7(1): e30136.
  * Durante et al. (2017) Durante, D., Dunson, D. B., and Vogelstein, J. T. (2017).  “Nonparametric Bayes modeling of populations of networks.”  Journal of the American Statistical Association, 112(520): 1516–1530.
  * Durante et al. (2018) Durante, D., Dunson, D. B., et al. (2018).  “Bayesian inference and testing of group differences in brain networks.”  Bayesian Analysis, 13(1): 29–58.
  * Everitt (2012) Everitt, R. G. (2012).  “Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks.”  Journal of Computational and Graphical Statistics, 21(4): 940–960.
URL https://doi.org/10.1080/10618600.2012.687493

  * Fallani et al. (2017) Fallani, F. D. V., Latora, V., and Chavez, M. (2017).  “A Topological Criterion for Filtering Information in Complex Brain Networks.”  PLOS Computational Biology, 13(1): e1005305.
URL
https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005305

  * Frank and Strauss (1986) Frank, O. and Strauss, D. (1986).  “Markov graphs.”  Journal of the american Statistical association, 81(395): 832–842.
  * Fuster (2006) Fuster, J. M. (2006).  “The cognit: A network model of cortical representation.”  International Journal of Psychophysiology, 60(2): 125 – 132.
URL http://www.sciencedirect.com/science/article/pii/S0167876006000079

  * Geerligs et al. (2017) Geerligs, L., Tsvetanov, K. A., Cam-CAN, and Henson, R. N. (2017).  “Challenges in measuring individual differences in functional connectivity using fMRI: The case of healthy aging.”  Human Brain Mapping, 38(8): 4125–4156.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/hbm.23653

  * Gelfand et al. (1995) Gelfand, A. E., Sahu, S. K., and Carlin, B. P. (1995).  “Efficient Parametrisations for Normal Linear Mixed Models.”  Biometrika, 82(3): 479–488.
URL http://www.jstor.org/stable/2337527

  * Geman and Geman (1984) Geman, S. and Geman, D. (1984).  “Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images.”  IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6(6): 721–741.
  * Ginestet et al. (2017) Ginestet, C. E., Li, J., Balachandran, P., Rosenberg, S., and Kolaczyk, E. D. (2017).  “Hypothesis testing for network data in functional neuroimaging.”  Ann. Appl. Stat., 11(2): 725–750.
URL https://doi.org/10.1214/16-AOAS1015

  * Ginestet et al. (2011) Ginestet, C. E., Nichols, T. E., Bullmore, E. T., and Simmons, A. (2011).  “Brain Network Analysis: Separating Cost from Topology Using Cost-Integration.”  PLOS ONE, 6(7): 1–17.
URL https://doi.org/10.1371/journal.pone.0021570

  * Haario et al. (2001) Haario, H., Saksman, E., Tamminen, J., et al. (2001).  “An adaptive Metropolis algorithm.”  Bernoulli, 7(2): 223–242.
  * Hastings (1970) Hastings, W. K. (1970).  “Monte Carlo Sampling Methods Using Markov Chains and Their Applications.”  Biometrika, 57(1): 97–109.
URL http://www.jstor.org/stable/2334940

  * Hoff et al. (2002) Hoff, P. D., Raftery, A. E., and Handcock, M. S. (2002).  “Latent space approaches to social network analysis.”  Journal of the american Statistical association, 97(460): 1090–1098.
  * Holland et al. (1983) Holland, P. W., Laskey, K. B., and Leinhardt, S. (1983).  “Stochastic blockmodels: First steps.”  Social networks, 5(2): 109–137.
  * Hunter et al. (2008a) Hunter, D., Handcock, M., Butts, C., Goodreau, S., and Morris, M. (2008a).  “ergm: A Package to Fit, Simulate and Diagnose Exponential-Family Models for Networks.”  Journal of Statistical Software, Articles, 24(3): 1–29.
URL https://www.jstatsoft.org/v024/i03

  * Hunter et al. (2008b) Hunter, D. R., Goodreau, S. M., and Handcock, M. S. (2008b).  “Goodness of Fit of Social Network Models.”  Journal of the American Statistical Association, 103(481): 248–258.
URL https://doi.org/10.1198/016214507000000446

  * Kolaczyk (2009) Kolaczyk, E. D. (2009).  Statistical Analysis of Network Data: Methods and Models.  Springer Series in Statistics. New York: Springer-Verlag.
URL //www.springer.com/gb/book/9780387881454

  * Koskinen (2004) Koskinen, J. (2004).  “Bayesian analysis of exponential random graphs-estimation of parameters and model selection.”  Technical report, Stockholm University.
  * Koskinen et al. (2013) Koskinen, J. H., Robins, G. L., Wang, P., and Pattison, P. E. (2013).  “Bayesian analysis for partially observed network data, missing ties, attributes and actors.”  Social Networks, 35(4): 514–527.
URL http://www.sciencedirect.com/science/article/pii/S0378873313000671

  * Krivitsky (2012) Krivitsky, P. N. (2012).  “Exponential-family random graph models for valued networks.”  Electron. J. Statist., 6: 1100–1128.
URL http://dx.doi.org/10.1214/12-EJS696

  * Krivitsky and Handcock (2014) Krivitsky, P. N. and Handcock, M. S. (2014).  “Supplementary Material: A separable model for dynamic networks.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1): 29–46.
URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/rssb.12014

  * Lehmann et al. (2021) Lehmann, B., Henson, R., Geerligs, L., Cam-CAN, and White, S. (2021).  “Characterising group-level brain connectivity: A framework using Bayesian exponential random graph models.”  NeuroImage, 225: 117480.
URL https://www.sciencedirect.com/science/article/pii/S1053811920309654

  * Mukherjee et al. (2017) Mukherjee, S. S., Sarkar, P., and Lin, L. (2017).  “On Clustering Network-Valued Data.”  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, 7074–7084. Red Hook, NY, USA: Curran Associates Inc.
  * Murray et al. (2006) Murray, I., Ghahramani, Z., and MacKay, D. J. C. (2006).  “MCMC for doubly-intractable distributions.”  In Proceedings of the 22nd Annual Conference on Uncertainty in Artificial Intelligence (UAI-06), 359–366. AUAI Press.
  * Obando and Fallani (2017) Obando, C. and Fallani, F. D. V. (2017).  “A statistical model for brain networks inferred from large-scale electrophysiological signals.”  Journal of The Royal Society Interface, 14(128): 20160940.
URL http://rsif.royalsocietypublishing.org/content/14/128/20160940

  * Papaspiliopoulos et al. (2003) Papaspiliopoulos, O., Roberts, G. O., and Sköld, M. (2003).  “Non-centered parameterisations for hierarchical models and data augmentation.”  In Bernardo, J., Bayarri, M., Berger, J., Dawid, A., Heckerman, D., Smith, A., and West, M. (eds.), Bayesian Statistics 7: Proceedings of the Seventh Valencia International Meeting, volume 307. Oxford University Press, USA.
  * Papaspiliopoulos et al. (2007) — (2007).  “A General Framework for the Parametrization of Hierarchical Models.”  Statist. Sci., 22(1): 59–73.
URL https://doi.org/10.1214/088342307000000014

  * Park and Haran (2018) Park, J. and Haran, M. (2018).  “Bayesian inference in the presence of intractable normalizing functions.”  Journal of the American Statistical Association, 113(523): 1372–1390.
  * Roberts et al. (1997) Roberts, G. O., Gelman, A., Gilks, W. R., et al. (1997).  “Weak convergence and optimal scaling of random walk Metropolis algorithms.”  Annals of Applied probability, 7(1): 110–120.
  * Roberts and Rosenthal (2009) Roberts, G. O. and Rosenthal, J. S. (2009).  “Examples of adaptive MCMC.”  Journal of Computational and Graphical Statistics, 18(2): 349–367.
  * Roberts et al. (2001) Roberts, G. O., Rosenthal, J. S., et al. (2001).  “Optimal scaling for various Metropolis-Hastings algorithms.”  Statistical science, 16(4): 351–367.
  * Schweinberger and Handcock (2015) Schweinberger, M. and Handcock, M. S. (2015).  “Local dependence in random graph models: characterization, properties and statistical inference.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77(3): 647–676.
URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/rssb.12081

  * Schweinberger et al. (2020) Schweinberger, M., Krivitsky, P. N., Butts, C. T., and Stewart, J. R. (2020).  “Exponential-Family Models of Random Graphs: Inference in Finite, Super and Infinite Population Scenarios.”  Statist. Sci., 35(4): 627–662.
URL https://doi.org/10.1214/19-STS743

  * Shafto et al. (2014) Shafto, M. A., Tyler, L. K., Dixon, M., Taylor, J. R., Rowe, J. B., Cusack, R., Calder, A. J., Marslen-Wilson, W. D., Duncan, J., Dalgleish, T., Henson, R. N., Brayne, C., and Matthews, F. E. (2014).  “The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing.”  BMC Neurology, 14(1): 1–25.
URL http://dx.doi.org/10.1186/s12883-014-0204-1

  * Signorelli and Wit (2020) Signorelli, M. and Wit, E. C. (2020).  “Model-based clustering for populations of networks.”  Statistical Modelling, 20(1): 9–29.
  * Simpson et al. (2013) Simpson, S., Lyday, R., Hayasaka, S., Marsh, A., and Laurienti, P. (2013).  “A permutation testing framework to compare groups of brain networks.”  Frontiers in Computational Neuroscience, 7: 171.
URL https://www.frontiersin.org/article/10.3389/fncom.2013.00171

  * Simpson et al. (2011) Simpson, S. L., Hayasaka, S., and Laurienti, P. J. (2011).  “Exponential random graph modeling for complex brain networks.”  PloS one, 6(5): e20039.
  * Sinke et al. (2016) Sinke, M. R., Dijkhuizen, R. M., Caimo, A., Stam, C. J., and Otte, W. M. (2016).  “Bayesian exponential random graph modeling of whole-brain structural networks across lifespan.”  NeuroImage, 135(Supplement C): 79 – 91.
URL http://www.sciencedirect.com/science/article/pii/S1053811916301069

  * Slaughter and Koehly (2016) Slaughter, A. J. and Koehly, L. M. (2016).  “Multilevel models for social networks: hierarchical Bayesian approaches to exponential random graph modeling.”  Social networks, 44: 334–345.
  * Stillman et al. (2017) Stillman, P. E., Wilson, J. D., Denny, M. J., Desmarais, B. A., Bhamidi, S., Cranmer, S. J., and Lu, Z.-L. (2017).  “Statistical Modeling of the Default Mode Brain Network Reveals a Segregated Highway Structure.”  Scientific Reports, 7(1): 11694.
URL https://doi.org/10.1038/s41598-017-09896-6

  * Sweet et al. (2013) Sweet, T. M., Thomas, A. C., and Junker, B. W. (2013).  “Hierarchical network models for education research: Hierarchical latent space models.”  Journal of Educational and Behavioral Statistics, 38(3): 295–318.
  * Sweet et al. (2014) — (2014).  “Hierarchical mixed membership stochastic blockmodels for multiple networks and experimental interventions.”  In Airoldi, E. M., Blei, D., Erosheva, E. A., and Fienberg, S. E. (eds.), Handbook of Mixed Membership Models and Their Applications, 463–488. Chapman & Hall/CRC, 1st edition.
  * Tan and Friel (2020) Tan, L. S. and Friel, N. (2020).  “Bayesian variational inference for exponential random graph models.”  Journal of Computational and Graphical Statistics, 1–19.
  * Thiemichen et al. (2016) Thiemichen, S., Friel, N., Caimo, A., and Kauermann, G. (2016).  “Bayesian exponential random graph models with nodal random effects.”  Social Networks, 46: 11–28.
URL http://www.sciencedirect.com/science/article/pii/S0378873316000034

  * Tierney (1994) Tierney, L. (1994).  “Markov Chains for Exploring Posterior Distributions.”  Ann. Statist., 22(4): 1701–1728.
URL https://doi.org/10.1214/aos/1176325750

  * Tzourio-Mazoyer et al. (2002) Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., Mazoyer, B., and Joliot, M. (2002).  “Automated Anatomical Labeling of Activations in SPM Using a Macroscopic Anatomical Parcellation of the MNI MRI Single-Subject Brain.”  NeuroImage, 15(1): 273 – 289.
URL http://www.sciencedirect.com/science/article/pii/S1053811901909784

  * Wang and Atchadé (2014) Wang, J. and Atchadé, Y. F. (2014).  “Approximate Bayesian computation for exponential random graph models for large social networks.”  Communications in Statistics-Simulation and Computation, 43(2): 359–377.
  * Wang et al. (2013) Wang, P., Robins, G., Pattison, P., and Lazega, E. (2013).  “Exponential random graph models for multilevel networks.”  Social Networks, 35(1): 96–115.
  * Yin et al. (2020) Yin, F., Shen, W., and Butts, C. T. (2020).  “Finite Mixtures of ERGMs for Modeling Ensembles of Networks.”
  * Yu and Meng (2011) Yu, Y. and Meng, X.-L. (2011).  “To Center or Not to Center: That Is Not the Question—An Ancillarity–Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC Efficiency.”  Journal of Computational and Graphical Statistics, 20(3): 531–570.
URL https://doi.org/10.1198/jcgs.2011.203main