# Bayesian model-based clustering for multiple network data

Anastasia Mantziou1 Simón Lunagómez1 Robin Mitra2

( 1Mathematics and Statistics, Lancaster University
2School of Mathematics, Cardiff University )

###### Abstract

There is increasing appetite for analysing multiple network data. This is
different to analysing traditional data sets, where now each observation in
the data comprises a network. Recent technological advancements have allowed
the collection of this type of data in a range of different applications. This
has inspired researchers to develop statistical models that most accurately
describe the probabilistic mechanism that generates a network population and
use this to make inferences about the underlying structure of the network
data. Only a few studies developed to date consider the heterogeneity that can
exist in a network population. We propose a Mixture of Measurement Error
Models for identifying clusters of networks in a network population, with
respect to similarities detected in the connectivity patterns among the
networks’ nodes. Extensive simulation studies show our model performs well in
both clustering multiple network data and inferring the model parameters. We
further apply our model on two real-world multiple network data sets resulting
from the fields of Computing (Human Tracking Systems) and Neuroscience.

Keywords: Bayesian models, Clustering, Mixture models, Multiple network data,
Object data analysis

##  1 Introduction

Conventional statistical methods for modelling and analysing data, such as
standard regression models, assume each observation or outcome value is a
scalar. While this is appropriate for many applications there is sometimes a
need to handle more complex types of data. For example, each observation may
constitute a set of interconnected points where the distribution of the
connections between the points determines the observation’s properties. Such
data are typically referred to as networks in the literature, and modelling
them is fundamentally important in some applications. For example, suppose we
are interested in monitoring movements of subjects across a geographical area
via displays at fixed locations. Each display represents a point, or node, of
the network and when a subject moves from one display to another a connection,
or edge, between the two nodes is assumed. The pattern of movement across the
displays then characterises a network for that subject. Modelling such data
presents some unique challenges. For example, we may be interested in
distinguishing different patterns of behaviour among the different subjects.
Are some subjects more likely to visit certain displays than others? Do
subjects have different or unusual patterns of movements to the majority of
subjects? These are all important applied questions that present significant
challenges to address due to the complexity of dealing with network data.

The availability of multiple network data has risen substantially in recent
years, due to the advancement of technological means that record this type of
data (White et al., (1986), Fields and Song, (1989)). This has inspired many
researchers to develop statistical models that most accurately describe the
probabilistic mechanism that generates a network population. Specifically,
there are three different frameworks considered in the literature for
modelling multiple network data: 1) the latent space models, 2) the distance-
based models, and 3) the measurement error models.

For a single network observation, the fundamental idea behind the latent space
class of models is that the occurrence of an edge between two nodes depends on
the positions of the nodes in a latent space (Hoff et al., (2002), Young and
Scheinerman, (2007)). Recent studies (Gollini and Murphy, (2016), Durante et
al., (2017), Wang et al., (2019), Nielsen and Witten, (2018), Arroyo et al.,
(2019)) on modelling multiple network data have extended this idea to build
models for populations of networks with aligned vertex sets, assuming that the
nodes lie in a common, unobserved subspace.

Another approach to modelling multiple network data is the utilisation of
distance metrics that measure similarities among networks with respect to
global or local characteristics of the networks (Donnat and Holmes, (2018)).
Under this framework, researchers rely on the notion of an average network
that represents a network population, with respect to a specified distance
metric (Lunagómez et al., (2019), Kolaczyk et al., (2017), Ginestet et al.,
(2017)).

The third class of models, measurement error models, account for the erroneous
nature of the networks. A fundamental source of noise found in network data
originates from the various measurement tools used for the construction of
networks, i.e. the processes used to measure an interaction (edge) between two
objects (nodes). Researchers focusing on the statistical analysis of networks
as single observations have developed methods to incorporate the uncertainty
of falsely observing edges or non-edges in a network. Such studies involve
predicting network topologies accounting for the falsely non-observed edges
(Jiang et al., (2011)), estimating the adjacency matrix from a set of noisy
entries (Chatterjee et al., (2015)), classifying nodes of networks with
errorful edges (Priebe et al., (2015)), developing a regression model for
networks assuming that the observed network is a perturbed version of a true
unobserved network (Le and Li, (2020)) and performing Bayesian inference on
the network’s structure utilising information from measurements (Young et al.,
(2020)). Another group of studies focuses on the propagation of the error to
network summary statistics (Balachandran et al., (2017), Chang et al.,
(2020)), and to estimators of average causal effects under network
interference when the error arises from a measurement process used to
construct the network (Li et al., (2021)). Le et al., (2018), Newman, (2018)
and Peixoto, (2018) develop this idea to model multiple network observations
in order to infer the probabilistic mechanism that generates the network
population. Specifically they assume the networks are noisy realisations of a
true unobserved network.

Despite the growing research interest on modelling multiple network data, only
few studies developed to date consider the heterogeneity that can exist in a
network population. Notably, Mukherjee et al., (2017) were the first to
consider the problem of clustering multiple network data. They assume two
different scenarios, (a) the networks in the population share the same set of
nodes, and (b) the networks in the population do not share the same set of
nodes. Our paper focuses on Scenario (a). In Scenario (a), the authors obtain
a mixture model of graphons and implement a spectral clustering algorithm to
infer the membership allocation of each network observation.

An application driven study on clustering multiple network data is introduced
by Diquigiovanni and Scarpa, (2019), who aim to cluster a population of
networks where each network observation represents the playing style of a
football team at a specific match. The clustering approach seen in this study
involves the specification of an ad hoc measure of similarity between
networks, and the implementation of an agglomerative method for clustering the
networks according to their similarities.

To the best of our knowledge, the third and last study that examines the
problem of clustering network populations is that of Signorelli and Wit,
(2020). In this study, the authors deal with the problem of clustering using a
mixture model whose components can be any statistical network model, under the
restriction that it can be specified as a Generalised Linear Model (GLM). For
estimating the parameters of their model, they implement an Expectation
Maximisation (EM) algorithm for a predefined number of clusters. To determine
the network model for their mixture, the authors propose the initial use of a
mixture of saturated network models, to reveal information about the structure
of the data at hand. The saturated network model assumes that each edge in
each network in the population is generated with some unique, unconstrained
probability.

Another group of studies that accounts for the heterogeneity in a set of
network observations are the studies that perform the task of network
classification. Some of these studies consider either specific network summary
measures (Prasad et al., (2015)), or vectorise only the important entries
(edges) of the adjacency matrix (Richiardi et al., (2011), Zhang et al.,
(2012)) to classify networks. Thus, they ignore the overall networks’
structure. In contrast to these studies, Relión et al., (2019) perform
prediction of the class membership of networks using a linear classifier with
the adjacency matrices of the networks as predictors. Their approach accounts
for the networks’ structure by using a penalty to select important nodes and
edges.

These contributions provide interesting approaches for identifying variations
between network data, but there are some key limitations associated with
these:

  * •

In the studies of Mukherjee et al., (2017), Diquigiovanni and Scarpa, (2019)
and Relión et al., (2019), the methods proposed are non model-based. Mukherjee
et al., (2017) and Diquigiovanni and Scarpa, (2019) propose algorithms that
detect underlying network clusters in the data and Relión et al., (2019)
predict class membership of the networks. In all three studies the groups of
networks identified cannot be interpreted using a parametric representation.
Interpretability of the different groups of networks in a population is
crucial in many applications in order to infer group specific properties and
differences.

  * •

While Signorelli and Wit, (2020) provide a model-based approach for clustering
multiple network data, the mixture components must conform to rigid modelling
assumptions. This means that only specific characteristics of the networks can
be inferred depending on what these model assumptions allow. It would be ideal
to have a framework that is flexible enough to incorporate different modelling
assumptions as deemed appropriate to application allowing the most
scientifically relevant inferences to be made.

  * •

In addition, Signorelli and Wit, (2020) propose to initially obtain a mixture
of saturated network models, thus resulting in an overly complex model with a
large number of parameters to estimate. This can substantially increase the
computational time needed for the EM algorithm to converge as well as
increasing the potential for non-convergence due to having to explore a very
high dimensional parameter space.

  * •

The supervised approach of Relión et al., (2019) requires a training data set
to predict the class. The class labels of the networks in the training data
set must be pre-specified, which can be restrictive for some network
applications for which we do not have a priori information about the networks.

To address these limitations, in this paper we propose a mixture model for
identifying clusters of networks in a network population, with respect to
similarities detected in the connectivity patterns of the networks’ nodes. We
consider the case when the networks in the population share the same set of
$n$ nodes, and each network could belong to one of, a predefined number of,
$C$ clusters. Inspired by the approach of Le et al., (2018), we adopt a
measurement error formulation, assuming networks lying within each of the $C$
clusters are noisy realisations of a true underlying network representative.
The attractive feature of this approach is that it decouples the statistical
model for the network data from the underlying cluster specific network
properties. We are thus able to provide a flexible model-based approach for
detecting clusters of networks in a network population, as well as interpret
these clusters with respect to our model parameterisation. Our framework is
also flexible enough to incorporate, and thus exploit, any underlying
assumptions about the structure of the networks within the clusters that are
of scientific interest or otherwise supported by the data.

Le et al., (2018) develop a model for multiple network observations assuming
that noisy network-valued observations arise from a true underlying adjacency
matrix. The inferential framework built in their study consists of two steps.
First, they use a Spectral Clustering algorithm to infer the community
structure formed by the nodes of the true underlying network, and second they
implement an EM algorithm to estimate the model parameters. An evident
limitation of their inferential framework is that their algorithm does not
simultaneously update the parameters of their model for the network data and
the parameters characterising the underlying network structure, as this would
require the development of new techniques. In addition, the assumption of a
sole true underlying adjacency matrix is quite restrictive, especially for a
large sample of networks where a degree of heterogeneity is expected.

We adopt a Bayesian modelling approach which provides some unique advantages
over previous approaches. In particular by utilising Markov Chain Monte Carlo
(MCMC) methods, we are able to simultaneously infer the cluster membership of
the networks, together with model parameters characterising the distribution
of the networks within each cluster as well as those that characterise the
structure of the underlying cluster specific network representatives. To best
of our knowledge, there is no coherent framework in the literature that
permits this type of complete inference from the network data.

Our framework is flexible enough to answer a diverse range of applied
questions with respect to the heterogeneity in a network population. These
include being able to detect clusters of networks as well as inferring key
different features between clusters through comparisons between the underlying
representatives. In addition, when interest is in identifying observations
that do not follow the distribution of the majority of the network data, the
framework can also be formulated to detect outlying network observations. The
approach will be illustrated through two different applied examples, one
involving monitoring movement of people across a University Campus and another
measuring individuals’ connectivity patterns across different regions of the
brain.

The remainder of this article is organised as follows. In Section 2, we
describe the applied examples that motivated the development of the methods.
In Section 3, we provide background to our modelling framework. In Section 4,
we develop the Bayesian formulation of a mixture of measurement error models
for network data, along with the MCMC scheme to make inferences. In Section 5,
we present simulation studies to assess the performance of our method for
various network sizes and sample sizes. In Section 6, we analyse two different
multiple network data examples to illustrate the broad applicability of our
methods. Lastly, in Section 7, we give some concluding remarks.

##  2 Two motivating examples

In this section we introduce two different applied examples that have
triggered research questions which we aim to answer in our study. While both
applications address very different areas, a common feature in both data sets
is the heterogeneity of the network data in the sample.

###  2.1 Data on movements of subjects across a University Campus

The first example comes from data collected on movement of people across
Lancaster University Campus in the UK. The study was performed by members of
the Computing and Communications department at the university (Shaw et al.,
(2018)).

A series of fixed displays are located across the campus. Individuals taking
part in the study installed a Tacita mobile application on their phone, and
whenever they pass one of these displays this application registers their
presence at that location. The application can also serve as a mean of
communication between a display and a viewer, in order for the viewers to be
able to see content relevant to their interests. Specifically, the viewer can
request what to see on the screen of the display, but also the display can
detect when a user is in its proximity in order to show content aligned with
their interests. Thus, the application records the consecutive displays
visited by the users, along with the time visited, and the type of content
shown.

Consequently, the Tacita data set serves as an example of multiple network
observations, as the movements of each individual can be represented by a
network. Thus, for each individual we can obtain a network where nodes
represent displays, and edges represent the movements of the user among the
displays. In this example each network corresponds to the aggregated movements
of an individual across the campus during different times of the day,
resulting in 120 network observations. The following questions arise:

  * •

Can we detect different patterns among the users’ movements?

  * •

Can we cluster the users according to their movements?

  * •

How informative can the clustering be for the users in our data?

###  2.2 Data measuring connectivity patterns across different regions of the
brain

The second example is a multiple network data set arising from the field of
Neuroscience. In this data example, connectivity patterns across different
regions of the brain were measured for 30 healthy individuals at resting
state. For each individual a series of 10 measurements were taken using
diffusion magnetic resonance imaging (dMRI). The measurements are represented
as networks with the nodes corresponding to fixed regions of the brain, and
edges denoting the connections recorded among those regions. Specifically, the
network data consist of 200 nodes (regions of the brain) according to the
CC200 atlas (Craddock et al., (2012)).

This data set has been also discussed in the study of Zuo et al., (2014),
Arroyo et al., (2019) and Lunagómez et al., (2019), with the latter two
studies analysing the data from a modelling perspective. Specifically, Arroyo
et al., (2019) investigate the ability of their method to identify differences
among individuals with respect to communities formed from the networks’ nodes,
while Lunagómez et al., (2019) assume unimodality of the probabilistic
mechanism that generates the network population and infer a representative
network for the population of individuals, according to a pre-specified
distance metric. None of these approaches seek to determine and interpret
clusters of networks.

Our goal in this application is to explore possible heterogeneity amongst the
networks. Questions include:

  * •

Can we identify clusters of individuals with respect to similarities found in
their connectivity patterns?

  * •

Can we interpret those clusters with respect to the postulated model
parameterisation.

  * •

Can we determine whether there are any individuals with uncommon brain
connectivity patterns?

These applied research questions have motivated the proposed mixture of
network measurement error models, described in detail in Section 6.

##  3 Background

A network can be represented as a graph $\mathcal{G}=(V,E)$, where
$V=\\{1,\ldots,n\\}$ represents the set of $n$ nodes and $E$ represents the
set of observed edges in $\mathcal{G}$, with $E\in\mathcal{E}_{n}$ and
$\mathcal{E}_{n}=\\{(i,j)|i,j\in V\\}$. A common mathematical network
representation is an $n\times n$ adjacency matrix $A_{\mathcal{G}}$, with the
$A_{\mathcal{G}}(i,j)$ entry of the matrix denoting the state of the $(i,j)$
edge. The $(i,j)^{th}$ element of the adjacency matrix for a graph with binary
edges is,

$A_{\mathcal{G}}(i,j)=\begin{cases}1,\text{ if an edge occurs between nodes i
and j,}\\\ 0,\text{ otherwise}.\end{cases}$

The adjacency matrix of an undirected graph with no self-loops is symmetric
with $A_{\mathcal{G}}(i,j)=A_{\mathcal{G}}(j,i)$ and $A_{\mathcal{G}}(i,i)=0$,
for $i,j\in\\{1,\ldots,n\\}$. By $\mathcal{G}_{1},\ldots,\mathcal{G}_{N}$ we
represent a population of $N$ graphs, with corresponding adjacency matrices
$A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}}$. In this study, we assume
that the networks in the population $\mathcal{G}_{1},\ldots,\mathcal{G}_{N}$
are undirected with no self-loops, and share the same set of $n$ nodes. We
note here that our methods can be equally applied to populations of directed
graphs.

When modelling a population of network observations it is often natural to
assume that the networks have been subjected to some noise or measurement
error during their construction. For example, in the application investigating
the movement of subjects across campus, it is possible that the Tacita mobile
application fails to register a subject at a display occasionally, while also
sometimes incorrectly registering a subject at a display, particularly when
displays are located fairly close together. This results in network data that
might have edges missing as well as edges recorded that should not be present.

Under this measurement error hypothesis, the researcher assumes that the
observed network data correspond to noisy realisations of a true underlying
network, which leads to recording some erroneous edges, due to the existence
of an underlying measurement error process. Le et al., (2018) were the first
to introduce this approach for modelling populations of networks, which has
inspired the proposed modelling approach.

Le et al., (2018) assume that the information contained in the network
population can be summarised by a representative network $\mathcal{G}^{*}$,
and a measurement error process that does not allow us to accurately observe
the representative network. Specifically, the authors assume a false positive
probability $P_{ij}$ of observing an edge between nodes $i,j$ in the $k^{th}$
network observation $\mathcal{G}_{k}$, given that there is no edge between the
same two nodes in the representative network $\mathcal{G}^{*}$; and
respectively, a false negative probability $Q_{ij}$ of not observing an edge
for the nodes $i,j$ in the $k^{th}$ network observation $\mathcal{G}_{k}$,
while there is an edge for the same two nodes in the representative network
$\mathcal{G}^{*}$. Thus, the entries of the matrices $P,\text{ }Q$, are the
false positive/negative probabilities of seeing/not seeing an edge
respectively between two nodes in the data.

The mathematical formulation of the above set-up is the following. Let
$A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}$ denote the adjacency matrices
for the network population, and $A_{\mathcal{G}^{*}}$ the adjacency matrix of
the representative network. The false positive and false negative
probabilities $P_{ij},Q_{ij}$ can be described as follows,

$\displaystyle\text{if }A_{\mathcal{G}^{*}}(i,j)=1,\text{ then
}A_{\mathcal{G}_{k}}(i,j)=\begin{cases}1,\quad\text{with prob }1-Q_{ij}\\\
0,\quad\text{with prob }Q_{ij}\end{cases};$ $\displaystyle\text{if
}A_{\mathcal{G}^{*}}(i,j)=0,\text{ then
}A_{\mathcal{G}_{k}}(i,j)=\begin{cases}1,\quad\text{with prob }P_{ij}\\\
0,\quad\text{with prob }1-P_{ij}\end{cases}.$ (1)

From (1) it follows that the probability of the occurrence or non-occurrence
of an edge between nodes $i,j$ in the $k^{th}$ network observation is,

$P(A_{\mathcal{G}_{k}}(i,j)|A_{\mathcal{G}^{*}}(i,j)=1,P_{ij},Q_{ij})=(1-Q_{ij})^{A_{\mathcal{G}_{k}}(i,j)}\cdot
Q_{ij}^{1-A_{\mathcal{G}_{k}}(i,j)},\text{ if }A_{\mathcal{G}^{*}}(i,j)=1;$
$P(A_{\mathcal{G}_{k}}(i,j)|A_{\mathcal{G}^{*}}(i,j)=0,P_{ij},Q_{ij})=P_{ij}^{A_{\mathcal{G}_{k}}(i,j)}\cdot(1-P_{ij})^{1-A_{\mathcal{G}_{k}}(i,j)},\text{
if }A_{\mathcal{G}^{*}}(i,j)=0.$

Le et al., (2018) treat the adjacency matrix of the representative network
$A_{\mathcal{G}^{*}}$ as a latent variable, while the false positive and false
negative probabilities $P_{ij}$ and $Q_{ij}$ are model parameters. Thus, the
likelihood of the representative network $A_{\mathcal{G}^{*}}$ given the
network data $A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}$, as seen in Le
et al., (2018), is

$\displaystyle\mathcal{L}(A_{\mathcal{G}^{*}};A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}})=\prod_{k=1}^{N}\prod_{(i,j):i<j}[(1-Q_{ij})^{A_{\mathcal{G}_{k}}(i,j)}\cdot
Q_{ij}^{1-A_{\mathcal{G}_{k}}(i,j)}\cdot
W_{ij}]^{A_{\mathcal{G}^{*}}(i,j)}\cdot$
$\displaystyle[P_{ij}^{A_{\mathcal{G}_{k}}(i,j)}\cdot(1-P_{ij})^{1-A_{\mathcal{G}_{k}}(i,j)}\cdot(1-W_{ij})]^{1-A_{\mathcal{G}^{*}}(i,j)},$

where $W_{ij}=\mathbb{E}A_{\mathcal{G}^{*}}(i,j)$ represents the probability
of observing an edge between nodes $i,j$ in $A_{\mathcal{G}^{*}}$.

Le et al., (2018) further assume that the nodes of the underlying true network
form communities that can be described by a Stochastic Block Model (SBM).
Under the SBM assumption, each node of the true network belongs to an
unobserved block $k\in\\{1,\cdots,K\\}$, and the probability of observing an
edge between two nodes depends on their block membership denoted by
$\\{b_{i}\\}_{i=1}^{n}$, with $b_{i}\in\\{1,\cdots,K\\}$. The probability of
observing an edge between nodes $(i,j)\text{ with }b_{i}=k,\text{ }b_{j}=l$ is
represented by $\theta_{kl}$.

The inference of the model parameters and the latent variable is conducted in
two stages. First, a Spectral Clustering algorithm is applied to reveal the
underlying block membership of the representative’s nodes, and second, an EM
algorithm is implemented to estimate the model parameters. While this
formulation has appealing features it would be ideal to have a coherent
modelling framework that can jointly infer block membership of the
representative’s nodes together with the parameters characterising the
distribution of the network data. In addition, using an EM algorithm to
estimate model parameters means that measures of uncertainty such as standard
errors rely on asymptotic approximations that may not be valid in many
applications, particularly when involving small samples sizes.

In the next section we propose a mixture of measurement error models inspired
by Le et al., (2018) for clustering heterogeneous network data. We adopt a
Bayesian framework that allows us to jointly infer the parameters of the
measurement error model as well as those characterising the underlying network
representatives corresponding to each cluster. In addition, the Bayesian
formulation is flexible enough to accommodate diverse modelling assumptions
for the network representatives.

##  4 A Mixture of Measurement Error models

In this section, we detail the formulation and implementation of the mixture
of measurement error models. We first describe the Bayesian formulation of the
measurement error model when there is only one cluster. We then extend this to
multiple clusters. Following this we describe how posterior samples can be
obtained using MCMC. Finally we describe a special case of this formulation
that can correspond to detecting outlying networks in the data.

###  4.1 Model formulation

To begin with we assume underlying the network data there is a latent
representative network with adjacency matrix denoted by $A_{\mathcal{G}^{*}}$.
In addition, we assume that the probability of observing a false positive or
false negative edge between two nodes in the network data is independent of
the pair of nodes considered. Thus, the false positive probability, $p$, and
false negative probability, $q$, can be viewed as scalars. Under this
specification, the probability mass function of the edge state between nodes
$(i,j)$ in network observation $k$, given $A_{\mathcal{G}^{*}}(i,j)$, $p$,
$q$, is

$\displaystyle
P(A_{\mathcal{G}_{k}}(i,j)|A_{\mathcal{G}^{*}}(i,j),p,q)=[(1-q)^{A_{\mathcal{G}_{k}}(i,j)}\cdot
q^{1-A_{\mathcal{G}_{k}}(i,j)}]^{A_{\mathcal{G}^{*}}(i,j)}\cdot$
$\displaystyle[p^{A_{\mathcal{G}_{k}}(i,j)}\cdot(1-p)^{1-A_{\mathcal{G}_{k}}(i,j)}]^{1-A_{\mathcal{G}^{*}}(i,j)}.$

Hence, the conditional probability of observing a network population
$A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}$ given $A_{\mathcal{G}^{*}}$,
$p$, $q$ is,

$\displaystyle
P(A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}}|p,q,A_{\mathcal{G}^{*}})=\prod_{k=1}^{N}\prod_{(i,j):i<j}P(A_{\mathcal{G}_{k}}(i,j)|A_{\mathcal{G}^{*}}(i,j),p,q)=$
$\displaystyle\prod_{k=1}^{N}\prod_{(i,j):i<j}((1-q)^{A_{\mathcal{G}_{k}}(i,j)}\cdot
q^{1-A_{\mathcal{G}_{k}}(i,j)})^{A_{\mathcal{G}^{*}}(i,j)}\cdot(p^{A_{\mathcal{G}_{k}}(i,j)}\cdot(1-p)^{1-A_{\mathcal{G}_{k}}(i,j)})^{1-A_{\mathcal{G}^{*}}(i,j)}.$

An advantage of the measurement error formulation is that the model
specification for the representative network $A_{\mathcal{G}^{*}}$ can vary
depending on the type of information the analyst wants to capture for the data
at hand. As previously discussed, Le et al., (2018) assume a SBM for the
network representative. For illustration we also assume an SBM structure for
the representative $A_{\mathcal{G}^{*}}$, but note that this can be easily
modified, e.g. reduced to a simpler model such as the Erdös-Rényi if supported
by the data. The SBM for the representative can be represented hierarchically
in the following way,

$\displaystyle
A_{\mathcal{G}^{*}}(i,j)|\boldsymbol{\theta},\boldsymbol{b}\sim\text{Bernoulli}(\theta_{b_{i}b_{j}});$
$\displaystyle\theta_{kl}\sim\text{Beta}(\epsilon_{kl},\zeta_{kl});$
$\displaystyle\boldsymbol{b|w}\sim\text{Multinomial}(\boldsymbol{w});$

where $w_{k}$ represents the probability of a node to belong to block
$k\in\\{1,\cdots,K\\}$. For the probability vector
$\boldsymbol{w}=\\{w_{1},\cdots,w_{K}\\}$ we assume a symmetric Dirichlet
prior distribution with hyperparameter $\boldsymbol{\chi}$. Common choices for
the hyperparameter vector $\boldsymbol{\chi}$ are setting all elements to 0.5
or 1.

Thus the hierarchical structure of the model is,

$\displaystyle\prod_{k=1}^{N}\prod_{(i,j):i<j}P(A_{\mathcal{G}_{k}}(i,j)|A_{\mathcal{G}^{*}}(i,j),p,q)P(A_{\mathcal{G}^{*}(i,j)}|\boldsymbol{\theta},\boldsymbol{b})$

where

$P(A_{\mathcal{G}^{*}}(i,j)|\boldsymbol{\theta},\boldsymbol{b})=\theta_{b_{i}b_{j}}^{A_{\mathcal{G}^{*}}(i,j)}(1-\theta_{b_{i}b_{j}})^{1-A_{\mathcal{G}^{*}}(i,j)}.$

We further specify a Beta prior distribution for both the false positive $p$
and false negative $q$ probabilities,

$p\sim\text{Beta}(\alpha_{0},\beta_{0})\text{,
}q\sim\text{Beta}(\gamma_{0},\delta_{0}),$

which facilitates posterior computations. A common choice sets the Beta prior
hyperparameters to 0.5, corresponding to the Jeffreys prior.

###  4.2 Mixture of measurement error models

We further extend the measurement error model to a mixture of measurement
error models, with a predefined number of mixture components, $C$, in order to
provide a model-based approach for identifying clusters of networks in a
network population $A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}$. We assume
each cluster $c$ of networks is described by a unique network representative
$A_{\mathcal{G}^{*}_{c}}$, a false positive probability $p_{c}$, and a false
negative probability $q_{c}$, where $c\in\\{1,\ldots,C\\}$. In this section,
we present the Bayesian framework for this mixture of measurement error
models. Each cluster-specific representative network is characterised by an
SBM, and the block structure of each representative is allowed to vary.

Let $\boldsymbol{z}=(z_{1},\cdots,z_{N})\in\\{1,\cdots,C\\}$ be the latent
variables representing the cluster membership of the network data
$A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}$. Then the conditional
probability of the data given $\boldsymbol{z}$ takes the form

$\displaystyle
P(A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}}|\\{p_{c},q_{c},A_{\mathcal{G}^{*}_{c}}\\}_{c=1}^{C},z_{1},\cdots,z_{N})=$
$\displaystyle\prod_{k=1}^{N}\Big{(}\prod_{(i,j):i<j}\Big{(}(1-q_{z_{k}})^{A_{\mathcal{G}_{k}}(i,j)}q_{z_{k}}^{1-A_{\mathcal{G}_{k}}(i,j)}\Big{)}^{A_{\mathcal{G}^{*}_{z_{k}}}(i,j)}\cdot\Big{(}p_{z_{k}}^{A_{\mathcal{G}_{k}}(i,j)}(1-p_{z_{k}})^{1-A_{\mathcal{G}_{k}}(i,j)}\Big{)}^{1-A_{\mathcal{G}^{*}_{z_{k}}}(i,j)}\Big{)}.$

The log of the conditional probability of the data up to a normalising
constant is

$\displaystyle\sum_{k=1}^{N}\Bigg{(}\sum_{(i,j):i<j}\bigg{(}A_{\mathcal{G}^{*}_{z_{k}}}(i,j)\Big{(}A_{\mathcal{G}_{k}}(i,j)\log\frac{(1-q_{z_{k}})(1-p_{z_{k}})}{q_{z_{k}}p_{z_{k}}}+\log\frac{q_{z_{k}}}{1-p_{z_{k}}}\Big{)}+$
$\displaystyle+\Big{(}A_{\mathcal{G}_{k}}(i,j)\log\frac{p_{z_{k}}}{1-p_{z_{k}}}+\log(1-p_{z_{k}})\Big{)}\bigg{)}\Bigg{)}.$

We assume that the cluster labels $z_{1},\cdots,z_{N}$ follow a Multinomial
distribution with parameter $\boldsymbol{\tau}=(\tau_{1},\cdots,\tau_{C})$,
where $\tau_{c}$ represents the probability that a network observation belongs
to cluster $c$, and $\sum_{c=1}^{C}\tau_{c}=1$. We assume a symmetric
Dirichlet prior distribution for the vector of probabilities
$\boldsymbol{\tau}$ which has the advantage of being conditionally conjugate
with the distribution for $\boldsymbol{z}$. As commented previously common
choices set the Dirichlet hyperparameters all to 0.5 or to 1.

###  4.3 MCMC scheme for mixture model

With the modelling framework described above we are able to draw samples from
the joint posterior distribution of the parameters using MCMC. The joint
posterior distribution is known up to a normalising constant, specifically

$\displaystyle\begin{gathered}P(\boldsymbol{A_{\mathcal{G}^{*}}},\boldsymbol{p},\boldsymbol{q},\boldsymbol{z},\boldsymbol{\tau},\boldsymbol{w},\boldsymbol{b},\boldsymbol{\theta}|A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}})\\\
\propto
P(A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}}|\boldsymbol{A_{\mathcal{G}^{*}}},\boldsymbol{p},\boldsymbol{q},\boldsymbol{z})\cdot
P(\boldsymbol{A_{\mathcal{G}^{*}}}|\boldsymbol{w},\boldsymbol{b},\boldsymbol{\theta})\cdot
P(\boldsymbol{p}|\boldsymbol{\alpha_{0}},\boldsymbol{\beta_{0}})\\\ \cdot
P(\boldsymbol{q}|\boldsymbol{\gamma_{0}},\boldsymbol{\delta_{0}})\cdot
P(\boldsymbol{z}|\boldsymbol{\tau})\cdot
P(\boldsymbol{\tau}|\boldsymbol{\psi})\cdot
P(\boldsymbol{\theta}|\boldsymbol{\epsilon},\boldsymbol{\zeta})\cdot
P(\boldsymbol{b}|\boldsymbol{w})\cdot
P(\boldsymbol{w}|\boldsymbol{\chi}),\end{gathered}$ (5)

where
$P(A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}}|\boldsymbol{A_{\mathcal{G}^{*}}},\boldsymbol{p},\boldsymbol{q},\boldsymbol{z})$
is the conditional probability of the network data,
$P(\boldsymbol{A_{\mathcal{G}^{*}}}|\boldsymbol{w},\boldsymbol{b},\boldsymbol{\theta})$
is the conditional probability of the latent variable
$\boldsymbol{A_{\mathcal{G}^{*}}}$ and the rest of the components of the right
hand side of the expression are the prior distributions for the model
parameters, as defined in Sections 4.1 and 4.2.

To obtain posterior samples from the joint posterior in (2) we note that many
of the full conditional distributions of the unknown quantities
(parameters/latent data) are available in closed form, and when these are not
available these can be approximated using Metropolis-Hastings. As a result we
obtain posterior inferences through a component wise MCMC sampler, also known
as a Metropolis-Hastings-within-Gibbs sampler. This closely follows a Gibbs
sampler, where all parameters and latent data are updated from their full
conditional distributions except for the full conditionals of
$\\{A_{\mathcal{G}^{*}_{c}},p_{c},q_{c}\\}_{c=1}^{C}$, which are approximated
using Metropolis-Hastings proposal distributions.

In the Metropolis-Hastings step, we use a mixture of kernels for updating the
parameters of the measurement error model
$\\{A_{\mathcal{G}^{*}_{c}},p_{c},q_{c}\\}_{c=1}^{C}$, in analogy to the MCMC
scheme seen in Lunagómez et al., (2019). Specifically in every iteration of
the MCMC we update the adjacency matrix of the network representative of
cluster $c$, $A_{\mathcal{G}^{*}_{c}}$, using either of the following two
proposals with some fixed probability:

  * (I)

We perturb the edges of the current network representative
$A_{\mathcal{G}^{*}_{c}}^{(curr)}$ of cluster $c$ in the following way:

$A_{\mathcal{G}^{*}_{c}}^{(prop)}(i,j)=\begin{cases}1-A_{\mathcal{G}^{*}_{c}}^{(curr)}(i,j),\text{
with probability }\omega\\\ A_{\mathcal{G}^{*}_{c}}^{(curr)}(i,j),\text{ with
probability }1-\omega\end{cases}.$

  * (II)

We propose a new network representative $A_{\mathcal{G}^{*}_{c}}^{(prop)}$ for
cluster $c$ drawing each edge of the proposed representative
$A_{\mathcal{G}^{*}_{c}}^{(prop)}(i,j)$ independently from a Bernoulli
distribution with parameter
$\frac{1}{N}\sum_{k=1}^{N}A_{\mathcal{G}_{k}}(i,j)$.

For case (I), we accept the proposed network representative
$A_{\mathcal{G}^{*}_{c}}^{(prop)}$ with probability

$\min\Bigg{\\{}1,\frac{P(A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}|A_{\mathcal{G}^{*}_{c}}^{(prop)},p_{c}^{(curr)},q_{c}^{(curr)},\boldsymbol{z}^{(curr)})P(A_{\mathcal{G}^{*}_{c}}^{(prop)}|\boldsymbol{b_{c}}^{(curr)},\boldsymbol{\theta_{c}}^{(curr)})}{P(A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}|A_{\mathcal{G}^{*}_{c}}^{(curr)},p_{c}^{(curr)},q_{c}^{(curr)},\boldsymbol{z}^{(curr)})P(A_{\mathcal{G}^{*}_{c}}^{(curr)}|\boldsymbol{b_{c}}^{(curr)},\boldsymbol{\theta_{c}}^{(curr)})}\Bigg{\\}},$
(6)

where
$P(A_{\mathcal{G}^{*}_{c}}^{(\cdot)}|\boldsymbol{b_{c}},\boldsymbol{\theta_{c}})$
is the SBM assumed for the representative defined in Section 4.1. The proposal
distribution under case (I) is symmetric, and so it cancels out from the
Metropolis ratio in expression (3).

Under case (II), we accept the proposed network representative
$A_{\mathcal{G}^{*}_{c}}^{(prop)}$ with probability

$\displaystyle\min\Bigg{\\{}1,\frac{P(A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}|A_{\mathcal{G}^{*}_{c}}^{(prop)},p_{c}^{(curr)},q_{c}^{(curr)},\boldsymbol{z}^{(curr)})P(A_{\mathcal{G}^{*}_{c}}^{(prop)}|\boldsymbol{b_{c}}^{(curr)},\boldsymbol{\theta_{c}}^{(curr)})}{P(A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}|A_{\mathcal{G}^{*}_{c}}^{(curr)},p_{c}^{(curr)},q_{c}^{(curr)},\boldsymbol{z}^{(curr)})P(A_{\mathcal{G}^{*}_{c}}^{(curr)}|\boldsymbol{b_{c}}^{(curr)},\boldsymbol{\theta_{c}}^{(curr)})}$
$\displaystyle\cdot\frac{Q(A_{\mathcal{G}^{*}_{c}}^{(curr)}|A_{\mathcal{G}^{*}_{c}}^{(prop)})}{Q(A_{\mathcal{G}^{*}_{c}}^{(prop)}|A_{\mathcal{G}^{*}_{c}}^{(curr)})}\Bigg{\\}},$
(7)

where $Q(A_{\mathcal{G}^{*}_{c}}^{(\cdot)}|A_{\mathcal{G}^{*}_{c}}^{(\cdot)})$
corresponds to the proposal distribution.

To update the false positive probability $p_{c}$ of cluster $c$, we use a
mixture of random walk proposals indexed by $l$ following Lunagómez et al.,
(2019)).

  * •

Draw $v\sim\text{Unif}(-u_{l},u_{l})$, for $0<u_{l}<0.5$.

  * •

Calculate the candidate proposal value $y=p_{c}^{(curr)}+v$.

  * •

Propose a new value for $p_{c}$ (constrained to lie in the interval (0,0.5)
for identifiability reasons) as follows,

$p_{c}^{(prop)}=\begin{cases}y,\quad\text{if }0<y<0.5;\\\ -y,\quad\text{if
}y<0;\\\ 1-y,\quad\text{if }y>0.5.\\\ \end{cases}$

The mixture is over $\\{u_{1},\ldots,u_{L}\\}$. Thus, we perturb the current
state of the false positive probability $p_{c}^{(curr)}$ using various sizes
of $u_{l}$, each imposing a less or more drastic change on $p_{c}^{(curr)}$.
We accept the proposed value $p_{c}^{(prop)}$ with probability

$\min\Bigg{\\{}1,\frac{P(A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}|A_{\mathcal{G}^{*}_{c}}^{(curr)},p_{c}^{(prop)},q_{c}^{(curr)},\boldsymbol{z}^{(curr)})P(p_{c}^{(prop)}|\alpha_{0,c},\beta_{0,c})}{P(A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}|A_{\mathcal{G}^{*}_{c}}^{(curr)},p_{c}^{(curr)},q_{c}^{(curr)},\boldsymbol{z}^{(curr)})P(p_{c}^{(curr)}|\alpha_{0,c},\beta_{0,c})}\Bigg{\\}},$
(8)

where $P(p_{c}^{(\cdot)}|\alpha_{0,c},\beta_{0,c})$ is a
Beta($\alpha_{0,c},\beta_{0,c}$) prior as in Section 4.1. The proposal
distribution for $p_{c}$ is symmetric, thus it does not appear in the
Metropolis ratio in expression (5). In exactly the same manner, we update the
false negative probability $q_{c}$, for $c\in\\{1,\ldots,C\\}$.

The rest of the parameters are updated via Gibbs samplers, by drawing values
from their full conditional posteriors.

The full conditional posterior for $\boldsymbol{\tau}$ is given by

$\displaystyle
P(\boldsymbol{\tau}|\boldsymbol{A_{\mathcal{G}^{*}}},\boldsymbol{p},\boldsymbol{q},\boldsymbol{z},\boldsymbol{w},\boldsymbol{b},\boldsymbol{\theta},A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}})=\text{Dirichlet}(\psi+\eta_{1},\ldots,\psi+\eta_{C}).$
(9)

where $\eta_{c}=\sum_{j=1}^{N}1_{c}(z_{j})$, $c=1,\ldots,C$, denotes the
number of networks that belong to cluster $c$.

We draw the latent cluster-membership $z_{k}$ for each network observation $k$
from a Multinomial distribution with unnormalised probabilities specified in
the following way:

$\displaystyle
P(z_{k}=c|\boldsymbol{\tau},\boldsymbol{A_{\mathcal{G}^{*}}},\boldsymbol{p},\boldsymbol{q},A_{\mathcal{G}_{k}})\propto
P(A_{\mathcal{G}_{k}}|z_{k}=c,\boldsymbol{A_{\mathcal{G}^{*}}},\boldsymbol{p},\boldsymbol{q})\cdot
P(z_{k}=c|\boldsymbol{\tau})$
$\displaystyle=\tau_{c}\cdot\prod_{(i,j):i<j}\Big{(}(1-q_{c})^{A_{\mathcal{G}_{k}}(i,j)}q_{c}^{1-A_{\mathcal{G}_{k}}(i,j)}\Big{)}^{A_{\mathcal{G}^{*}_{c}}(i,j)}\cdot\Big{(}p_{c}^{A_{\mathcal{G}_{k}}(i,j)}(1-p_{c})^{1-A_{\mathcal{G}_{k}}(i,j)}\Big{)}^{1-A_{\mathcal{G}^{*}_{c}}(i,j)}$
(10)

where
$P(A_{\mathcal{G}_{k}}|z_{k}=c,\boldsymbol{A_{\mathcal{G}^{*}}},\boldsymbol{p},\boldsymbol{q})$
is the probability we observe network $k$ given cluster membership $z_{k}=c$,
described by a measurement error model. The normalised probabilities are
obtained via Bayes Theorem.

The full conditional posterior for $\boldsymbol{w_{c}}$ is

$P(\boldsymbol{w_{c}}|A_{\mathcal{G}^{*}_{c}},p_{c},q_{c},\boldsymbol{z},\boldsymbol{\tau},\boldsymbol{b_{c}},\boldsymbol{\theta_{c}},A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}})=\text{Dirichlet}(\chi+h_{1},\ldots,\chi+h_{K}).$

where $h_{k}$ denotes the number of the nodes that belong to block k.

The full conditional posterior for the vector of the block-specific
probabilities of an edge occurrence for the network representative of cluster
$c$, $\boldsymbol{\theta_{c}}$, is

$\displaystyle
P(\boldsymbol{\theta_{c}}|A_{\mathcal{G}^{*}_{c}},p_{c},q_{c},\boldsymbol{z},\boldsymbol{\tau},\boldsymbol{b_{c}},\boldsymbol{w_{c}},A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}})=\text{Beta}(A_{\mathcal{G}^{*}_{c}}[kl]+\epsilon_{kl},\zeta_{kl}+n_{c,kl}-A_{\mathcal{G}^{*}_{c}}[kl]).$
(11)

where
$A_{\mathcal{G}^{*}_{c}}[kl]=\sum_{(i,j):b_{c,i}=k,b_{c,j}=l}A_{\mathcal{G}^{*}_{c}}(i,j)$
represents the sum of the entries for the pairs of nodes of the network
representative for cluster $c$ that have block membership $k,l$ respectively,
and $n_{c,kl}=\sum_{(i,j):i\neq j}\mathbb{I}(b_{c,i}=k,b_{c,j}=l)$ represents
the number of the pair of nodes of the representative for cluster $c$ that
have membership $k,l$ respectively.

Similarly to the formulation obtained for updating the latent cluster-
membership $\boldsymbol{z}$ of the network data, we obtain updates of the
latent block-membership $\boldsymbol{b_{c}}$ for the nodes of the network
representative of cluster $c$ from a Multinomial distribution with
unnormalised probabilities specified as follows:

$\displaystyle
P(b_{c,i}=k|A_{\mathcal{G}^{*}_{c}},p_{c},q_{c},\boldsymbol{z},\boldsymbol{\tau},\boldsymbol{\theta_{c}},\boldsymbol{w_{c}},A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}})\propto
P(A_{\mathcal{G}^{*}_{c}}|\boldsymbol{w_{c}},\boldsymbol{\theta_{c}},b_{c,i}=k)\cdot
P(b_{c,i}=k|\boldsymbol{w_{c}})$
$\displaystyle=w_{c,k}\cdot\prod_{j=1}^{n}\theta_{kb_{c,j}}^{A_{\mathcal{G}^{*}_{c}}(i,j)}(1-\theta_{kb_{c,j}})^{1-A_{\mathcal{G}^{*}_{c}}(i,j)}.$
(12)

where
$P(A_{\mathcal{G}^{*}_{c}}|\boldsymbol{w_{c}},\boldsymbol{\theta_{c}},b_{c,i}=k)$
is the probability of observing the representative of cluster $c$,
$A_{\mathcal{G}^{*}_{c}}$, described by an SBM, given its $i^{th}$ node
belongs to block $k$. Normalised probabilities are obtained by Bayes Theorem.

For the detailed derivation of the full conditional posterior distributions
refer to the Supplementary material Section 1. In addition, the MCMC algorithm
for clustering is sketched in the Supplementary material Algorithm 1.

###  4.4 Outlier network detection

In this section we present a modification of the mixture model presented in
Section 4.2 that allows us to explore the heterogeneity in a population of
networks with a different perspective. Notably, we modify our mixture model to
detect a cluster of outlier networks that are different to the majority of the
networks in the population. Under this formulation, we are able to address
additional applied research questions from the data at hand.

In contrast to the mixture model formulated for multiple cluster
representatives, the outlier cluster detection model assumes a single network
representative for the whole population of networks. Under this setup, we
assume that there are ultimately two clusters of networks formed within the
population of networks, one cluster being the majority cluster, while the
other cluster determining the outlier networks in the population. Thus, while
the false positive and false negative probabilities remain component specific
for each of the two clusters, the network representative is no longer a
component specific latent variable.

Similar to the mixture model formulated in Section 4.2, we now specify the
number of clusters to $C=2$, and
$\boldsymbol{z}=(z_{1},\cdots,z_{N})\in\\{1,2\\}$ denotes the latent cluster
membership of the network data
$A_{\mathcal{G}_{1}},\ldots,A_{\mathcal{G}_{N}}$. Under the assumption of a
single network representative, $A_{\mathcal{G}^{*}}$, the conditional
probability of the data given the latent variables, $\boldsymbol{z}$, takes
the form,

$\displaystyle
P(A_{\mathcal{G}_{1}},\cdots,A_{\mathcal{G}_{N}}|\\{p_{c},q_{c}\\}_{c=1}^{C},A_{\mathcal{G}^{*}},z_{1},\cdots,z_{N})=$
$\displaystyle\prod_{k=1}^{N}\Big{(}\prod_{(i,j):i<j}\Big{(}(1-q_{z_{k}})^{A_{\mathcal{G}_{k}}(i,j)}q_{z_{k}}^{1-A_{\mathcal{G}_{k}}(i,j)}\Big{)}^{A_{\mathcal{G}^{*}}(i,j)}\cdot\Big{(}p_{z_{k}}^{A_{\mathcal{G}_{k}}(i,j)}(1-p_{z_{k}})^{1-A_{\mathcal{G}_{k}}(i,j)}\Big{)}^{1-A_{\mathcal{G}^{*}}(i,j)}\Big{)}.$

Again, the model specification for the representative network can vary
depending on the type of information we want to capture for the data at hand.
A common choice is to consider an SBM structure again for the representative.
Due to having now a single representative only, the SBM model parameters are
no more component specific to the cluster.

To sample from the joint posterior of this model, we develop a Metropolis-
Hastings-within-Gibbs MCMC scheme, as presented in Section 4.3. The full
conditional posterior distributions are obtained as seen in Section 4.3 with
the only difference that the parameters/latent variables characterising the
representative, $A_{\mathcal{G}^{*}}$, namely, the block-specific edge
probabilities, $\boldsymbol{\theta}$, the probability of a node to belong to a
block, $\boldsymbol{w}$, and the block membership of the nodes,
$\boldsymbol{b}$, are no longer component specific, i.e. not indexed by
cluster $c$.

##  5 Simulations

In this Section, we perform simulation studies to assess the performance of
our algorithm in inferring the model parameters/latent variables and
clustering network data. First, we explore the performance of our algorithm
for moderate network sizes and various noise levels and SBM models, and
second, we investigate the algorithm performance for various network and
sample sizes.

###  5.1 Moderate network sizes

In this simulation study we investigate the performance of our model in
inferring model parameters for network populations with a moderate size of
nodes in different scenarios. Specifically, we consider the case of networks
with $n=21$ nodes, and a population of $N=180$ networks. We assume $C=3$
clusters of networks in the population, and consider SBMs for each
representative network with $B=2$ blocks. We vary the model parameter values
in order to explore performance.

To simulate the network population, we first simulate the representative
network of each cluster. We generate representatives of the clusters under two
different SBM structures with the following parameters,

  * •

SBM 1: $(w_{1}^{(1)},w_{2}^{(1)})=(0.5,0.5)$,
$(w_{1}^{(2)},w_{2}^{(2)})=(0.5,0.5)$, $(w_{1}^{(3)},w_{2}^{(3)})=(0.5,0.5)$,
$(\theta_{11}^{(1)},\theta_{12}^{(1)},\theta_{22}^{(1)})=(0.8,0.2,0.8)$,
$(\theta_{11}^{(2)},\theta_{12}^{(2)},\theta_{22}^{(2)})=(0.8,0.2,0.8)$,
$(\theta_{11}^{(3)},\theta_{12}^{(3)},\theta_{22}^{(3)})=(0.8,0.2,0.8)$.

  * •

SBM 2: $(w_{1}^{(1)},w_{2}^{(1)})=(0.7,0.3)$,
$(w_{1}^{(2)},w_{2}^{(2)})=(0.5,0.5)$, $(w_{1}^{(3)},w_{2}^{(3)})=(0.3,0.7)$,
$(\theta_{11}^{(1)},\theta_{12}^{(1)},\theta_{22}^{(1)})=(0.7,0.05,0.8)$,
$(\theta_{11}^{(2)},\theta_{12}^{(2)},\theta_{22}^{(2)})=(0.7,0.05,0.8)$,
$(\theta_{11}^{(3)},\theta_{12}^{(3)},\theta_{22}^{(3)})=(0.7,0.05,0.8)$.

In Figure 1 we visualise the 21-node representatives for each of the $C=3$
clusters under each of the SBM structures (1 and 2) described above.

Figure 1: Top: Network representatives for clusters $c=1,$ $2$ and $3$
respectively (from left to right), under SBM structure 1. Bottom: Network
representatives for clusters $c=1,$ $2$ and $3$ respectively (from left to
right), under SBM structure 2.

Next, we generate a population of 180 networks by perturbing the edges of each
representative through a measurement error process. Specifically, we generate
edges for 60 networks in each of the $C=3$ clusters, depending on the
existence or non-existence of an edge in the representative network of the
corresponding cluster $c$, given a false positive $p_{c}$ and false negative
$q_{c}$ probability. The simulation regimes considered for $p_{c},q_{c}$ and
the SBM parameters described above, are presented in Table 1.

Sim | $\boldsymbol{(p_{1},p_{2},p_{3})}$ | $\boldsymbol{(q_{1},q_{2},q_{3})}$ | $\boldsymbol{(w_{1}^{(1)},w_{2}^{(1)})}$ | $\boldsymbol{(w_{1}^{(2)},w_{2}^{(2)})}$ | $\boldsymbol{(w_{1}^{(3)},w_{2}^{(3)})}$
---|---|---|---|---|---
1 | (0.1,0.1,0.1) | (0.2,0.2,0.2) | (0.5,0.5) | (0.5,0.5) | (0.5,0.5)
2 | (0.1,0.1,0.1) | (0.2,0.2,0.2) | (0.7,0.3) | (0.5,0.5) | (0.3,0.7)
3 | (0.1,0.1,0.1) | (0.3,0.3,0.3) | (0.5,0.5) | (0.5,0.5) | (0.5,0.5)
4 | (0.1,0.1,0.1) | (0.3,0.3,0.3) | (0.7,0.3) | (0.5,0.5) | (0.3,0.7)
5 | (0.2,0.2,0.2) | (0.1,0.1,0.1) | (0.5,0.5) | (0.5,0.5) | (0.5,0.5)
6 | (0.2,0.2,0.2) | (0.1,0.1,0.1) | (0.7,0.3) | (0.5,0.5) | (0.3,0.7)
7 | (0.2,0.2,0.2) | (0.3,0.3,0.3) | (0.5,0.5) | (0.5,0.5) | (0.5,0.5)
8 | (0.2,0.2,0.2) | (0.3,0.3,0.3) | (0.7,0.3) | (0.5,0.5) | (0.3,0.7)
9 | (0.3,0.3,0.3) | (0.1,0.1,0.1) | (0.5,0.5) | (0.5,0.5) | (0.5,0.5)
10 | (0.3,0.3,0.3) | (0.1,0.1,0.1) | (0.7,0.3) | (0.5,0.5) | (0.3,0.7)
11 | (0.3,0.3,0.3) | (0.2,0.2,0.2) | (0.5,0.5) | (0.5,0.5) | (0.5,0.5)
12 | (0.3,0.3,0.3) | (0.2,0.2,0.2) | (0.7,0.3) | (0.5,0.5) | (0.3,0.7)

Sim | $\boldsymbol{(\theta_{11}^{(1)},\theta_{12}^{(1)},\theta_{22}^{(1)})}$ | $\boldsymbol{(\theta_{11}^{(2)},\theta_{12}^{(2)},\theta_{22}^{(2)})}$ | $\boldsymbol{(\theta_{11}^{(3)},\theta_{12}^{(3)},\theta_{22}^{(3)})}$
---|---|---|---
1 | (0.8,0.2,0.8) | (0.8,0.2,0.8) | (0.8,0.2,0.8)
2 | (0.7,0.05,0.8) | (0.7,0.05,0.8) | (0.7,0.05,0.8)
3 | (0.8,0.2,0.8) | (0.8,0.2,0.8) | (0.8,0.2,0.8)
4 | (0.7,0.05,0.8) | (0.7,0.05,0.8) | (0.7,0.05,0.8)
5 | (0.8,0.2,0.8) | (0.8,0.2,0.8) | (0.8,0.2,0.8)
6 | (0.7,0.05,0.8) | (0.7,0.05,0.8) | (0.7,0.05,0.8)
7 | (0.8,0.2,0.8) | (0.8,0.2,0.8) | (0.8,0.2,0.8)
8 | (0.7,0.05,0.8) | (0.7,0.05,0.8) | (0.7,0.05,0.8)
9 | (0.8,0.2,0.8) | (0.8,0.2,0.8) | (0.8,0.2,0.8)
10 | (0.7,0.05,0.8) | (0.7,0.05,0.8) | (0.7,0.05,0.8)
11 | (0.8,0.2,0.8) | (0.8,0.2,0.8) | (0.8,0.2,0.8)
12 | (0.7,0.05,0.8) | (0.7,0.05,0.8) | (0.7,0.05,0.8)

Table 1: Simulation regimes for 21-node networks and $C=3$ clusters.

For each simulation regime, we run our MCMC for 500,000 iterations with a
burn-in of 150,000. We assess the performance of our model in inferring the
model parameters by obtaining the posterior means and credible intervals for
these parameters. For simulation regime 4, the posterior means obtained for
the false positive probabilities of each cluster $(p_{1},p_{2},p_{3})$ are
$(0.11,0.10,0.10)$, while the credible intervals of each false positive
probability are $(0.10,0.11)$, $(0.09,0.10)$ and $(0.09,0.11)$ respectively.
In addition, for the same simulation regime, the posterior means of the false
negative probabilities of each cluster $(q_{1},q_{2},q_{3})$ are
$(0.30,0.31,0.30)$, while the credible intervals obtained are $(0.29,0.31)$,
$(0.30,0.32)$ and $(0.29,0.32)$ respectively. For simulation regime 7, that
involves the highest levels of noise, the posterior means obtained for the
false positive probabilities of each cluster $(p_{1},p_{2},p_{3})$ are
$(0.20,0.20,0.19)$ and the credible intervals for each cluster specific false
positive probability are $(0.19,0.21)$,$(0.19,0.21)$ and $(0.18,0.20)$.
Respectively, for the same simulation regime the posterior mean of the false
negative probabilities $(q_{1},q_{2},q_{3})$ are $(0.29,0.29,0.30)$, and the
credible interval for each false negative probability respectively are
$(0.28,0.30)$, $(0.28,0.31)$ and $(0.28,0.31)$. We observe that the posterior
means obtained for the false positive and false negative probabilities are
very close to their true values. This finding also holds for the rest of the
simulation regimes. For the complete results for all the regimes, refer to the
Supplementary Material, Section 2.1, Tables 1-13.

In order to investigate the performance of our algorithm in identifying the
true representatives, we obtain the Hamming distance between the posterior
representative samples and the true representatives, after a burn-in of
150,000 and a lag of 50, leaving 7,000 posterior samples. We calculate the
proportion of times that the distance is less or equal to 1, 5 and 10
respectively, similar to the summaries obtained in Lunagómez et al., (2019).
For each simulation regime, we observe that all the posterior representative
samples drawn for each cluster have a Hamming distance from the true
representative less than or equal to 1, 5, and 10, 100% of the time, as
presented in the Supplementary material, Section 2.1, Table 14. This result
suggests that the true representatives are almost fully identified from our
algorithm.

In addition, we assess the effectiveness of our algorithm in identifying the
cluster membership of the networks using the clustering entropy and purity
indices. We obtain 7,000 posterior draws (after a burn-in of 150,000 and lag
of 50) for the cluster membership $z$, calculate the clustering entropy and
clustering purity with respect to the true membership of the networks, and
calculate their mean for each simulation regime. We note that a clustering
entropy value of 0 and clustering purity value of 1 indicate a perfect cluster
allocation of the networks. The simulation results indicate the true cluster
membership of the networks is fully recovered by our MCMC algorithm, with mean
entropy 0 and mean purity 1 for each simulation regime. These results are in
the Supplementary material, Section 2.1, Table 15.

We further investigate the performance of the model for various noise levels
in the network population. Specifically, we generate network populations by
perturbing the edges of the representatives of SBM structure 1, illustrated in
Figure 1, with varying sizes of the false positive and negative probabilities.
Specifically we consider the simulation regimes presented in Tables 3 and 3.

Sim | $\boldsymbol{(p_{1},p_{2},p_{3})}$ | $\boldsymbol{(q_{1},q_{2},q_{3})}$
---|---|---
1 | (0.01,0.01,0.01) | (0.1,0.1,0.1)
2 | (0.05,0.05,0.05) | (0.1,0.1,0.1)
3 | (0.1,0.1,0.1) | (0.1,0.1,0.1)
4 | (0.15,0.15,0.15) | (0.1,0.1,0.1)
5 | (0.2,0.2,0.2) | (0.1,0.1,0.1)
6 | (0.25,0.25,0.25) | (0.1,0.1,0.1)
7 | (0.3,0.3,0.3) | (0.1,0.1,0.1)
8 | (0.35,0.35,0.35) | (0.1,0.1,0.1)
9 | (0.4,0.4,0.4) | (0.1,0.1,0.1)
10 | (0.45,0.45,0.45) | (0.1,0.1,0.1)
Table 2: Simulation regimes for varying sizes of false positive probabilities
$p_{c}$ and fixed false negative probabilities $q_{c}$, for $c\in\\{1,2,3\\}$
and 21-node networks.

Sim | $\boldsymbol{(p_{1},p_{2},p_{3})}$ | $\boldsymbol{(q_{1},q_{2},q_{3})}$
---|---|---
1 | (0.1,0.1,0.1) | (0.01,0.01,0.01)
2 | (0.1,0.1,0.1) | (0.05,0.05,0.05)
3 | (0.1,0.1,0.1) | (0.1,0.1,0.1)
4 | (0.1,0.1,0.1) | (0.15,0.15,0.15)
5 | (0.1,0.1,0.1) | (0.2,0.2,0.2)
6 | (0.1,0.1,0.1) | (0.25,0.25,0.25)
7 | (0.1,0.1,0.1) | (0.3,0.3,0.3)
8 | (0.1,0.1,0.1) | (0.35,0.35,0.35)
9 | (0.1,0.1,0.1) | (0.4,0.4,0.4)
10 | (0.1,0.1,0.1) | (0.45,0.45,0.45)
Table 3: Simulation regimes for varying sizes of false negative probabilities
$q_{c}$ and fixed false positive probabilities $p_{c}$, for $c\in\\{1,2,3\\}$
and 21-node networks.

Figure 2: Posterior means for false positive probabilities $p_{c}$ for
$c\in\\{1,2,3\\}$ (y axis), plotted against the true values of $p_{c}$ (x
axis). The red dotted line corresponds to the $y=x$.

Figure 3: Posterior means for false negative probabilities $q_{c}$ for
$c\in\\{1,2,3\\}$ (y axis), plotted against the true values of $q_{c}$ (x
axis). The red dotted line corresponds to the $y=x$.

For each regime presented in Tables 3 and 3, we generate a network population
and run the MCMC for 500,000 iterations with a burn-in of 150,000 iterations.
We obtain the posterior means for the varying false positive probability
(Table 3) and the varying false negative probability (Table 3) for each
cluster $c\in\\{1,2,3\\}$, and plot them against their true values as seen in
Figures 3 and 3 respectively.

We observe the posterior means lie mostly close to the $y=x$ line, indicating
that our model performs well in inferring the true false positive and false
negative probabilities, even for high noise levels. However, in Figure 3, for
the highest noise value of $p_{c}=0.45$ for $c\in\\{1,2,3\\}$, the posterior
mean of the false negative probability of cluster 3 $p_{3}$ is equal to 0.2,
substantially different to its true value. These results suggest that our
model performs well in most cases, even for high noise levels, but we must be
cautious when making inferences for network populations with great variability
in their structure.

###  5.2 Varying sizes of networks and network populations

We now explore how well our model infers the parameters with respect to
various network sizes and sample sizes. We keep $C=3$ clusters and $B=2$
blocks. We consider four different network sizes of 25, 50, 75 and 100 nodes,
and simulate populations of 45, 90, 135, 180, 225, 270 and 315 networks, for
each network size respectively.

To simulate the network populations, we first generate the network
representatives of each cluster assuming each follows an SBM. We specify the
parameters of the SBMs so that the expected degree of the resulting network
representatives is preserved for all network sizes. The representatives
obtained for cluster $c=1$ for the different network sizes are visualised in
Figures 7, 7, 7 and 7. In the Supplementary material, Section 2.2, Figures
1-4, we illustrate the rest of the representatives for $c=2,3$, for the
network sizes considered. The resulting populations are generated by
perturbing the edges of each network representative, for each network size
considered, with a false positive $p_{c}$ and false negative $q_{c}$
probability fixed at 0.08, for $c\in\\{1,2,3\\}$.

Figure 4: 25-node representative of cluster labelled 1.

Figure 5: 50-node representative of cluster labelled 1.

Figure 6: 75-node representative of cluster labelled 1.

Figure 7: 100-node representative of cluster labelled 1.

Figure 8: Left: Absolute error (y axis) for model parameters $\boldsymbol{p}$,
$\boldsymbol{q}$ and $\boldsymbol{\tau}$, for 25-node networks and varying
population sizes (x axis). Right: Absolute error (y axis) for model parameters
$\boldsymbol{p}$, $\boldsymbol{q}$ and $\boldsymbol{\tau}$, for 50-node
networks and varying population sizes (x axis).

Figure 9: Left: Absolute error (y axis) for model parameters $\boldsymbol{p}$,
$\boldsymbol{q}$ and $\boldsymbol{\tau}$, for 75-node networks and varying
population sizes (x axis). Right: Absolute error (y axis) for model parameters
$\boldsymbol{p}$, $\boldsymbol{q}$ and $\boldsymbol{\tau}$, for 100-node
networks and varying population sizes (x axis).

For each simulated data set we run the MCMC for 500,000 iterations with a
burn-in of 150,000. We demonstrate the performance of our model by obtaining
the absolute error of the model parameters for each simulation regime, as seen
in Figures 9 and 9. Specifically, the plots demonstrate the variability of the
absolute error (y axis) for various sample sizes (x axis). The absolute error
is the absolute value of the difference between the posterior means obtained
after burn-in, and the true value of the parameter. The different coloured
lines and dots correspond to the different clusters considered. The multiple
dots of the same color correspond to 5 replications of the MCMC, each on a
different randomly generated data set. The lines connect the absolute errors
averaged across the replications.

The plots indicate that the sample size affects the performance of our model
for the network sizes considered. We observe that as the number of nodes
increase, there is an increase in the required number of networks to preserve
the same level of accuracy in estimation. Nevertheless, even for large
networks and small sample sizes, it is encouraging to see the posterior means
obtained are not far away from their true values.

We further illustrate the performance of our model in identifying the true
representatives of each cluster for the various sample and network sizes, in
the Supplementary Material, Section 2.2, Figure 5.

##  6 Data examples

In this Section, we present the application of our mixture model on the two
real-world multiple network data sets presented in Section 2.

###  6.1 Movement patterns across campus

As introduced in Section 2.1, Tacita is a mobile phone application that
records the displays visited by users, along with the time visited and the
type of content shown on the display. One way to represent the data collected
from the Tacita application is through a network, where nodes correspond to
displays, and edges correspond to movements of users among the displays.
Consequently, we obtain a multiple network data set where each network
observation corresponds to a user’s movements across displays. The final data
sample consists of 120 undirected and unweighted network observations that
share the same set of 37 nodes corresponding to the displays across campus.
Names of the displays are presented in the Supplementary Material, Section
3.1, Table 16. As our mixture model requires the pre-specification of the
number of clusters $C$, we conduct exploratory data analysis (EDA) to explore
the potential number of clusters in the population. We considered various
network distance metrics, and for each metric we obtain a distance matrix that
contains the pairwise distances of the networks in the population. We obtain
the Multi Dimensional Scaling (MDS) plot for each distance matrix obtained.
The MDS algorithm maps objects in a 2-d space, respecting their pairwise
distances. The MDS plots obtained from the EDA are presented in the
Supplementary Material, Section 3.1, Figures 6 and 7.

The EDA results suggest that $C=3$ clusters of networks is reasonable. To
meaningfully initialise the networks’ cluster membership, we combine the
results from four different distance metrics; the Hamming, the Jaccard, the
$l_{2}$, and the wavelets metrics. For a descriptive review on the distance
metrics refer to Donnat and Holmes, (2018). Specifically, we use a k-memoid
algorithm using the R package kmed Budiaji, (2019) to determine four different
cluster memberships, corresponding to the four different metrics considered,
and determine the final cluster membership initialisation using majority vote,
i.e. by determining which networks are consistently allocated to one of the
three clusters among the four memberships obtained. We initialise the
representative of each cluster by generating its edges using independent
Bernoulli draws. The probability with which we draw an edge between two
specific nodes of the representative corresponds to the proportion of times
that we see that edge in the network data of the corresponding cluster. We
then initialise the nodes’ block membership of the representatives using SBM
estimates from the R package blockmodels INRA and Leger, (2015) that suggest
the presence of two underlying blocks. We note here that a simpler network
model, namely the Erdös-Rényi model, could also be applied to describe the
representative networks. We run our MCMC for 500,000 iterations with a burn-in
of 100,000.

In Figure 10, the left and middle plots present the proportion of times that
the nodes of the representative networks of clusters 2 and 3 respectively,
belong in each of the two blocks specified. The results for the representative
of cluster $c=1$ are presented in the Supplementary material, Section 3.1,
Figure 10. We note here, that a block structure is not identified for the
representative of cluster $c=1$. From Figure 10, we observe a similar block
membership is revealed for the representatives of clusters 2 and 3. However,
we notice differences in the block allocation of some nodes between the two
representatives, namely nodes labelled 1, 3, 5, 7, 10, 11, 12 and 13. In
addition, for the representative of cluster $c=3$, the nodes are more clearly
allocated to each of the two blocks, as seen from the proportions.

In Figure 10, the right plot corresponds to the proportion of times that an
individual is allocated to the clusters. We see most individuals are clearly
allocated to one of the three clusters. The proportion of the individuals
allocated in each cluster further supports the assumption of $C=3$ clusters.
Specifically, we notice clusters 2 and 3 contain only 17% and 18% of the
observations respectively, while cluster 1 contains approximately 65%.
Assuming a greater number of clusters could potentially lead to clusters
having only few network observations, making the inference of a network
representative for those clusters impractical.

Figure 10: Left: Proportion of times (y axis) that each node (x axis) of the
representative network of cluster labelled 2 is allocated in Blocks 1 or 2,
after a burn-in of 100,000 iterations. Middle: Proportion of times (y axis)
that each node (x axis) of the representative network of cluster labelled 3 is
allocated in Blocks 1 or 2, after a burn-in of 100,000 iterations. Right:
Proportion of times (y axis) an individual’s network (x axis) is allocated to
each of the 3 clusters, after a burn-in of 100,000 iterations.

In addition, in Figure 11 we obtain the network representatives of each
cluster, with node labels corresponding to the numbering seen in Table 16
given in the Supplementary material, Section 3.1, and layout similar to the
true location of the displays on campus. The two different colours of the
nodes denote the block membership inferred for each representative. In the
Supplementary Material, Section 3.1, Figures 8 and 9, we present trace plots
of the false negative and false positive probabilities for each of the three
clusters.

Figure 11: Left: Posterior mode for the representative of cluster 1, with
posterior mass 100%, for the last 100,000 iterations. Middle: Posterior mode
for the representative of cluster 2, with posterior mass 37%, for the last
100,000 iterations. Right: Posterior mode for the representative of cluster 3,
with posterior mass 15%, for the last 100,000 iterations.

For cluster $c=1$, we observe that the representative concentrating the whole
posterior mass is sparse having only few edges and no SBM structure.
Specifically, we note that the edges of the representative correspond to
movements of individuals among displays that are very close to each other
(e.g. edge between nodes 4 and 10 corresponding to displays both located in
the same building). In addition, most of the networks in the population are
allocated to this cluster with a very small false positive probability with
posterior mean 0.003 and very large false negative probability with posterior
mean 0.49. The small false positive probability indicates that the edges
observed in the network data are correctly recorded, while the high false
negative probability indicates a high possibility of edges in the network data
that we do not observe. This is a reasonable finding as we would anticipate
that the Tacita application might have missed movements of users among
displays due to WiFi connection issues. This remark is also justified by the
movements observed in the representative network that correspond to movements
among displays within the same building, suggesting that when a user is in a
building, the WiFi connection is preserved, and so the application can record
the movements of the user.

For cluster $c=2$, we observe that the posterior mode of the representative
network concentrates a relatively smaller posterior mass of 37%, while being
slightly denser compared to the representative of cluster 1. Moreover, an SBM
structure is revealed, with the mostly connected nodes belonging in the same
block. This representative reveals a specific movement pattern of the users
allocated in cluster 2, corresponding to the displays located at the central
part of the campus (nodes 5, 6, 16, 25, 26, 27), as well as displays located
at Infolab (nodes 1 and 2), and Furness College (4 and 10). However, there is
a smaller proportion of networks allocated to cluster 2 compared to cluster 1.
However, the small false positive (posterior mean 0.02) and high false
negative probability (posterior mean 0.49) again indicates that there might be
movements of users not recorded by the application.

Lastly, for cluster $c=3$ the posterior mode of the network representative
concentrates a smaller posterior mass compared to the other two
representatives equal to 15%. We notice similarities both in the block
structure and the connectivity patterns of the representatives of clusters 2
and 3. However, the representative of cluster 3 is notably denser, and some
new movements of individuals at displays labelled 3 and 13 (Faraday College
and County College) are discovered. We also note that for cluster 3, the
posterior means obtained for the false positive and false negative
probabilities are similar to cluster 2. In addition, clusters 2 and 3 have a
similar proportion of individuals. Overall, a common movement pattern is
discovered for individuals in clusters 2 and 3, indicating movements among
displays located in the central part of the campus.

###  6.2 Connectivity patterns in the brain

As described in Section 2.2, this example involves a population of 300
undirected networks corresponding to 10 brain-scans taken for 30 healthy
individuals via diffusion magnetic resonance imaging (dMRI). The nodes of the
networks correspond to regions of the brain, and edges denote connections
recorded among these regions. Specifically, the network data consist of 200
nodes according to the CC200 atlas (Craddock et al., (2012)). Our goal is to
identify a cluster of individuals with brain connectivity patterns that differ
compared to the majority of the networks in the population, and characterise
these with respect to a model parameterisation. We thus implement our outlier
cluster detection algorithm, as discussed in Section 4.4.

We initialise our algorithm similarly to the initialisation performed for the
Tacita application in Section 6.1. Hence, we determine an initial membership
of networks in two different clusters by implementing a k-memoid algorithm
using the R package kmed Budiaji, (2019). This is done for three distance
matrices that correspond to the Jaccard, the wavelets and the $l_{2}$ distance
metrics. We combine the results using majority vote, to obtain the initial
cluster membership of each network. We initialise the network representative
by generating its edges through independent Bernoulli draws, with
probabilities equal to the proportion of times the corresponding edge is
observed in the network data. We initialise the block membership of the
representative’s nodes by SBM estimation on the initial representative using
the R package blockmodels INRA and Leger, (2015). For the rest of the
parameters of the model we consider three different random initialisations,
and run the MCMC for each initialisation for 1,000,000 iterations.

In Figure 12 we present the trace plots of the false positive and false
negative probabilities for the majority cluster for all iterations under three
different initialisations. In the Supplementary Material, Section 3.2, Figure
11, we also present the trace plots of the false positive and false negative
probabilities for the outlier cluster. We observe that under the three
different initialisations the algorithm converges very quickly to the same
region. This is encouraging and suggests that a high posterior region has been
identified.

Figure 12: Left: Trace plot for false positive probability $p$ for majority
cluster for 1,000,000 iterations and three different initialisations. Right:
Trace plot for false negative probability $q$ for majority cluster for
1,000,000 iterations and three different initialisations.

We also compare the results from the three initialisations by obtaining the
posterior mode of the representative network for the last 50,000 iterations.
In Figure 13 (left) we obtain the posterior mode for the network
representative under the first initialisation (posterior mass of 0.08). The
colours of the nodes correspond to the block membership. In Figure 13
(middle), we present the not in common edges between the posterior modes of
the first and second initialisation to facilitate comparisons. The black edges
correspond to the edges present in the posterior mode of the first
initialisation and not present in the posterior mode of the second
initialisation, and the pink edges correspond to the edges present in the
posterior mode under the second initialisation and not present in the
posterior mode of the first initialisation. In Figure 13 (right), we similarly
present the not in common edges between the posterior modes of the first and
third initialisation. The posterior mode of the second initialisation has
posterior mass of 0.08, while the posterior mode of the third initialisation
has posterior mass of 0.16.

Figure 13: Left: Posterior mode of representative network from $1^{\text{st}}$
initialisation. Middle: Network with not in common edges between posterior
modes of representatives from $1^{\text{st}}$ and $2^{\text{nd}}$
initialisation. Right: Network with not in common edges between posterior
modes of representatives from $1^{\text{st}}$ and $3^{\text{rd}}$
initialisation. The nodes’ colours correspond to the block structure
identified under each initialisation.

There are three interesting findings with respect to the representative
inferred under the three different initialisations. First, there is only a
small proportion of edges not in common among the posterior modes of the three
different initialisations, considering the density of the graphs. Second, our
algorithm infers the same block structure for the three posterior modes of the
representative networks. Third, the posterior masses for the posterior mode
representatives are small, but is expected due to the high dimensional space
of the networks. In general the results are very encouraging given the size of
the networks considered.

Figure 14: Proportion of times that each of the 10 brain scans from the 30
individuals is allocated in cluster labelled 1 and/or 2.

We further observe a smaller posterior mean for the false negative probability
in the outlier cluster (0.32) compared to the majority cluster (0.37),
suggesting that the edges not observed in the network data of the outlier
cluster are more likely to be correct compared to the majority cluster. This
finding also suggests that the network data in the outlier cluster are sparser
compared to the majority cluster. Also, the small posterior means of the false
positive probability for both clusters (0.016 for the majority cluster and
0.025 for the outlier cluster) indicate that the edges observed in the network
data are likely correct.

In Figure 14, we present the results for the cluster membership $z$ of the
network observations for the first initialisation. Specifically, we calculate
the proportion of times that a network observation is allocated to the
majority or the outlier cluster, labelled by 1 or 2 respectively, from the
last 100,000 iterations. Each subfigure in Figure 14 shows the cluster
allocation of the 10 brain scans obtained for the same individual. We see that
our model mostly allocates scans of the same individual to the same cluster.
Thus, our model detects similarities among the brain scans of the same
individual, giving credence to our model clustering the networks sensibly.

This is a common finding with Arroyo et al., (2019) who performed semi-
supervised classification on the brain network population. However, our
approach is different to Arroyo et al., (2019) in two ways. First, we
implement an unsupervised method to infer underlying clusters of networks. We
only pre-define the number of cluster in the population. Second, the
interpretation of the results of our model-based clustering method compared to
Arroyo et al., (2019) differs significantly, as it reveals a cluster of
individuals whose brain connectivity patterns are different to a majority
group, and are interpreted through a parametric model.

##  7 Discussion

In this paper we introduced a mixture model for multiple network data that
allows us to identify clusters of networks in a population. To achieve this,
we formulated a mixture of measurement error models, and developed a Bayesian
framework that allows us to make inferences for all model parameters jointly.
This framework permits a diverse specification of the network model for the
representative networks in the population, determined according to the type of
information we want to exploit based on the data.

Through extensive simulations, we observed our method reliably inferred the
model parameters and cluster membership for moderate-sized networks, even for
regimes with high noise levels. This is an interesting result, as for high
noise levels there is great variability in the structure of the simulated
network population, making inference a challenging task. The results suggested
that our model can perform well for a range of real data applications.
Simulations also examined the model performance for large network and
population sizes. This had not been explored by Signorelli and Wit, (2020) who
also develop a model-based approach for clustering multiple network data. We
observed that the absolute errors of the posterior means for the parameters
were small, even for large networks and a relatively small sample size. This
suggests that our algorithm does not require a large number of network
observations to accurately infer the model parameters.

The clustering performed on the Tacita application revealed three different
movement patters of the users. This can be deduced by the network
representatives inferred for the clusters, which are primarily characterised
by their density. The cluster described by the sparser representative, reveals
movements among closely located displays on campus. The second and third
cluster identified, enclose denser networks in the population, and the
representative of each cluster reveals a specific movement pattern of the
individuals therein. As the majority of the networks are very sparse it was
encouraging to see that our model was able to separate out the denser networks
and further distinguish two different clusters among this subset of
individuals.

Analysis of the brain network data led to some interesting findings. First,
the results suggest we identified a high posterior region for the false
positive and false negative probabilities of each cluster. Second, a similar
posterior mode for the network representative has been inferred under three
different initialisations. These are especially encouraging given the high
dimensional space spanned by 200-node network representatives. Lunagómez et
al., (2019) who also obtain a network representative in terms of a Fréchet
mean for the same brain network population resorted to divide and conquer
methods to be able to make inferences, which was not required here.

Our model could potentially incorporate covariates at the node or edge level
to inform the inference. In some network applications, it is common to have
additional information about the nodes or the edges of the network. For
example, the Tacita mobile application also records the type of content shown
by the display at the time visited by the user. The incorporation of this
additional information could potentially lead to interesting additional
findings. It would be also interesting to perform a follow up analysis on data
recorded by the Tacita application after the emergence of the COVID-19
pandemic. This would allow us to investigate whether there is a change in the
movement patterns of users before and after the pandemic.

In the brain network example, although we assume one majority and one outlier
cluster, our framework can be more flexible. We could assume the presence of
two or more clusters of outlier networks that could describe different levels
of variability in the data. However, the data may not always support such
inferences.

An interesting result from both the Tacita and the brain networks application,
is the high false negative probabilities inferred for the clusters, attributed
to the networks’ sparsity. Network sparsity is a common issue in many real
world network applications. One way to deal with the sparsity would be to
consider shrinkage priors e.g. formulating the Horseshoe priors (Carvalho et
al., (2009)). Another way to account for network sparsity is to assume the
networks are partially observed. In the literature, partially observed
networks have been considered under two different perspectives. The coarsening
approach focuses on incorporating the coarsening mechanism, that allows us to
only partially observe the networks, in the model, and efficiently impute the
partially observed data thereafter (Heitjan and Rubin, (1991), Handcock and
Gile, (2010), Heitjan and Rubin, (1990), Kim and Hong, (2012)). Another
approach focuses on the missingness of certain edges and performs edge
prediction (Koskinen et al., (2013), Marchette and Hohman, (2015), Zhao et
al., (2017), Airoldi and Blocker, (2013)). Under the first approach, one way
we could incorporate the coarsening mechanism is through the assumption of a
sampling design, while under the second approach we could have a two-stage
method which would first involve performing link prediction, and second
performing inference. All the aforementioned approaches require significant
modifications of our model and present interesting avenues for future
research.

Finally, a natural extension of our model, is to allow inference for the
number of clusters as well. To achieve this, methods in the literature such as
the reversible jump MCMC (Richardson and Green, (1997)), and Dirichlet Process
(DP) mixture models (Neal, (2000)) could be adopted. Reversible jump MCMC
would be particularly challenging due to the MCMC moving between very
different dimensions to make inferences. DP mixture models would require a
different MCMC scheme but present an interesting possibility to consider.

Our paper contributes to the growing literature on clustering complex types of
data. Examples of such studies are the study of Lu et al., (2014) who focus on
identifying clusters of juggling cycles using Functional Principal Component
Analysis (FPCA) on diverse data objects, Song et al., (2007) who perform
model-based clustering on time-dependent gene expression data using Functional
Data Analysis, and Shen et al., (2013) who propose an agglomerative clustering
approach for clustering shape data according to their structure. In summary,
the flexibility of our modelling framework has been shown to address diverse
applied research questions with the potential to be widely applicable to many
fields.

## References

  * Airoldi and Blocker, (2013) Airoldi, E. M. and Blocker, A. W. (2013).  Estimating latent processes on a network from indirect measurements.  Journal of the American Statistical Association, 108(501):149–164.
  * Arroyo et al., (2019) Arroyo, J., Athreya, A., Cape, J., Chen, G., Priebe, C. E., and Vogelstein, J. T. (2019).  Inference for multiple heterogeneous networks with a common invariant subspace.  arXiv preprint arXiv:1906.10026.
  * Balachandran et al., (2017) Balachandran, P., Kolaczyk, E. D., and Viles, W. D. (2017).  On the propagation of low-rate measurement error to subgraph counts in large networks.  The Journal of Machine Learning Research, 18(1):2025–2057.
  * Budiaji, (2019) Budiaji, W. (2019).  kmed: Distance-Based K-Medoids.  R package version 0.3.0.
  * Carvalho et al., (2009) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2009).  Handling sparsity via the horseshoe.  In Artificial Intelligence and Statistics, pages 73–80. PMLR.
  * Chang et al., (2020) Chang, J., Kolaczyk, E. D., and Yao, Q. (2020).  Estimation of subgraph densities in noisy networks.  Journal of the American Statistical Association, pages 1–14.
  * Chatterjee et al., (2015) Chatterjee, S. et al. (2015).  Matrix estimation by universal singular value thresholding.  Annals of Statistics, 43(1):177–214.
  * Craddock et al., (2012) Craddock, R. C., James, G. A., Holtzheimer III, P. E., Hu, X. P., and Mayberg, H. S. (2012).  A whole brain fmri atlas generated via spatially constrained spectral clustering.  Human brain mapping, 33(8):1914–1928.
  * Diquigiovanni and Scarpa, (2019) Diquigiovanni, J. and Scarpa, B. (2019).  Analysis of association football playing styles: An innovative method to cluster networks.  Statistical Modelling, 19(1):28–54.
  * Donnat and Holmes, (2018) Donnat, C. and Holmes, S. (2018).  Tracking network dynamics: A survey of distances and similarity metrics.  arXiv preprint arXiv:1801.07351.
  * 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.
  * Fields and Song, (1989) Fields, S. and Song, O.-k. (1989).  A novel genetic system to detect protein–protein interactions.  Nature, 340(6230):245–246.
  * Ginestet et al., (2017) Ginestet, C. E., Li, J., Balachandran, P., Rosenberg, S., Kolaczyk, E. D., et al. (2017).  Hypothesis testing for network data in functional neuroimaging.  The Annals of Applied Statistics, 11(2):725–750.
  * Gollini and Murphy, (2016) Gollini, I. and Murphy, T. B. (2016).  Joint modeling of multiple network views.  Journal of Computational and Graphical Statistics, 25(1):246–265.
  * Handcock and Gile, (2010) Handcock, M. S. and Gile, K. J. (2010).  Modeling social networks from sampled data.  The Annals of Applied Statistics, 4(1):5.
  * Heitjan and Rubin, (1990) Heitjan, D. F. and Rubin, D. B. (1990).  Inference from coarse data via multiple imputation with application to age heaping.  Journal of the American Statistical Association, 85(410):304–314.
  * Heitjan and Rubin, (1991) Heitjan, D. F. and Rubin, D. B. (1991).  Ignorability and coarse data.  The annals of statistics, pages 2244–2253.
  * 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.
  * INRA and Leger, (2015) INRA and Leger, J.-B. (2015).  blockmodels: Latent and Stochastic Block Model Estimation by a ’V-EM’ Algorithm.  R package version 1.1.1.
  * Jiang et al., (2011) Jiang, X., Gold, D., and Kolaczyk, E. D. (2011).  Network-based auto-probit modeling for protein function prediction.  Biometrics, 67(3):958–966.
  * Kim and Hong, (2012) Kim, J. K. and Hong, M. (2012).  Imputation for statistical inference with coarse data.  Canadian Journal of Statistics, 40(3):604–618.
  * Kolaczyk et al., (2017) Kolaczyk, E., Lin, L., Rosenberg, S., and Walters, J. (2017).  Averages of unlabeled networks: Geometric characterization and asymptotic behavior.  arXiv preprint arXiv:1709.02793.
  * 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.
  * Le et al., (2018) Le, C. M., Levin, K., Levina, E., et al. (2018).  Estimating a network from multiple noisy realizations.  Electronic Journal of Statistics, 12(2):4697–4740.
  * Le and Li, (2020) Le, C. M. and Li, T. (2020).  Linear regression and its inference on noisy network-linked data.  arXiv preprint arXiv:2007.00803.
  * Li et al., (2021) Li, W., Sussman, D. L., and Kolaczyk, E. D. (2021).  Causal inference under network interference with noise.  arXiv preprint arXiv:2105.04518.
  * Lu et al., (2014) Lu, X., Marron, J., et al. (2014).  Analysis of juggling data: Object oriented data analysis of clustering in acceleration functions.  Electronic Journal of Statistics, 8(2):1842–1847.
  * Lunagómez et al., (2019) Lunagómez, S., Olhede, S. C., and Wolfe, P. J. (2019).  Modeling network populations via graph distances.  arXiv preprint arXiv:1904.07367.
  * Marchette and Hohman, (2015) Marchette, D. J. and Hohman, E. L. (2015).  Utilizing covariates in partially observed networks.  In 2015 18th International Conference on Information Fusion (Fusion), pages 166–172. IEEE.
  * Mukherjee et al., (2017) Mukherjee, S. S., Sarkar, P., and Lin, L. (2017).  On clustering network-valued data.  Advances in neural information processing systems.
  * Neal, (2000) Neal, R. M. (2000).  Markov chain sampling methods for dirichlet process mixture models.  Journal of computational and graphical statistics, 9(2):249–265.
  * Newman, (2018) Newman, M. E. (2018).  Estimating network structure from unreliable measurements.  Physical Review E, 98(6):062321.
  * Nielsen and Witten, (2018) Nielsen, A. M. and Witten, D. (2018).  The multiple random dot product graph model.  arXiv preprint arXiv:1811.12172.
  * Peixoto, (2018) Peixoto, T. P. (2018).  Reconstructing networks with unknown and heterogeneous errors.  Physical Review X, 8(4):041011.
  * Prasad et al., (2015) Prasad, G., Joshi, S. H., Nir, T. M., Toga, A. W., Thompson, P. M., (ADNI, A. D. N. I., et al. (2015).  Brain connectivity and novel network measures for alzheimer’s disease classification.  Neurobiology of aging, 36:S121–S131.
  * Priebe et al., (2015) Priebe, C. E., Sussman, D. L., Tang, M., and Vogelstein, J. T. (2015).  Statistical inference on errorfully observed graphs.  Journal of Computational and Graphical Statistics, 24(4):930–953.
  * Relión et al., (2019) Relión, J. D. A., Kessler, D., Levina, E., and Taylor, S. F. (2019).  Network classification with applications to brain connectomics.  The annals of applied statistics, 13(3):1648.
  * Richardson and Green, (1997) Richardson, S. and Green, P. J. (1997).  On bayesian analysis of mixtures with an unknown number of components (with discussion).  Journal of the Royal Statistical Society: series B (statistical methodology), 59(4):731–792.
  * Richiardi et al., (2011) Richiardi, J., Eryilmaz, H., Schwartz, S., Vuilleumier, P., and Van De Ville, D. (2011).  Decoding brain states from fmri connectivity graphs.  Neuroimage, 56(2):616–626.
  * Shaw et al., (2018) Shaw, P., Mikusz, M., Nurmi, P., and Davies, N. (2018).  Tacita: A privacy preserving public display personalisation service.  In Proceedings of the 2018 ACM International Joint Conference and 2018 International Symposium on Pervasive and Ubiquitous Computing and Wearable Computers, pages 448–451.
  * Shen et al., (2013) Shen, W., Wang, Y., Bai, X., Wang, H., and Latecki, L. J. (2013).  Shape clustering: Common structure discovery.  Pattern Recognition, 46(2):539–550.
  * Signorelli and Wit, (2020) Signorelli, M. and Wit, E. C. (2020).  Model-based clustering for populations of networks.  Statistical Modelling, 20(1):9–29.
  * Song et al., (2007) Song, J. J., Lee, H.-J., Morris, J. S., and Kang, S. (2007).  Clustering of time-course gene expression data using functional data analysis.  Computational biology and chemistry, 31(4):265–274.
  * Wang et al., (2019) Wang, S., Arroyo, J., Vogelstein, J. T., and Priebe, C. E. (2019).  Joint embedding of graphs.  IEEE transactions on pattern analysis and machine intelligence.
  * White et al., (1986) White, J. G., Southgate, E., Thomson, J. N., and Brenner, S. (1986).  The structure of the nervous system of the nematode caenorhabditis elegans.  Philos Trans R Soc Lond B Biol Sci, 314(1165):1–340.
  * Young et al., (2020) Young, J.-G., Cantwell, G. T., and Newman, M. (2020).  Bayesian inference of network structure from unreliable data.  Journal of Complex Networks, 8(6):cnaa046.
  * Young and Scheinerman, (2007) Young, S. J. and Scheinerman, E. R. (2007).  Random dot product graph models for social networks.  pages 138–149.
  * Zhang et al., (2012) Zhang, J., Cheng, W., Wang, Z., Zhang, Z., Lu, W., Lu, G., and Feng, J. (2012).  Pattern classification of large-scale functional brain networks: identification of informative neuroimaging markers for epilepsy.  PloS one, 7(5):e36733.
  * Zhao et al., (2017) Zhao, Y., Wu, Y.-J., Levina, E., and Zhu, J. (2017).  Link prediction for partially observed networks.  Journal of Computational and Graphical Statistics, 26(3):725–733.
  * Zuo et al., (2014) Zuo, X.-N., Anderson, J. S., Bellec, P., Birn, R. M., Biswal, B. B., Blautzik, J., Breitner, J. C., Buckner, R. L., Calhoun, V. D., Castellanos, F. X., et al. (2014).  An open science resource for establishing reliability and reproducibility in functional connectomics.  Scientific data, 1(1):1–13.