# Modelling Populations of Interaction Networks via Distance Metrics

George Bolt<EMAIL_ADDRESS>
STOR-i Centre for Doctoral Training
Lancaster University, Lancaster, UK, LA1 4YF Simón Lunagómez
<EMAIL_ADDRESS>
Department of Statistics
Instituto Tecnológico Autónomo de México (ITAM)
Río Hondo 1, Altavista, Álvaro Obregón, 01080 Ciudad de México, CDMX, Mexico
Christopher Nemeth<EMAIL_ADDRESS>
Department of Mathematics and Statistics
Lancaster University, Lancaster, UK, LA1 4YF

###### Abstract

Network data arises through observation of relational information between a
collection of entities. Recent work in the literature has independently
considered when (i) one observes a sample of networks, connectome data in
neuroscience being a ubiquitous example, and (ii) the units of observation
within a network are edges or paths, such as emails between people or a series
of page visits to a website by a user, often referred to as interaction
network data. The intersection of these two cases, however, is yet to be
considered. In this paper, we propose a new Bayesian modelling framework to
analyse such data. Given a practitioner-specified distance metric between
observations, we define families of models through location and scale
parameters, akin to a Gaussian distribution, with subsequent inference of
model parameters providing reasoned statistical summaries for this non-
standard data structure. To facilitate inference, we propose specialised
Markov chain Monte Carlo (MCMC) schemes capable of sampling from doubly-
intractable posterior distributions over discrete and multi-dimensional
parameter spaces. Through simulation studies we confirm the efficacy of our
methodology and inference scheme, whilst its application we illustrate via an
example analysis of a location-based social network (LSBN) data set.

Keywords: Interaction network, multiple networks, distance metrics, exchange
algorithm, involutive MCMC.

##  1 Introduction

Network data, defined to be information regarding relations amongst some
collection of entities, can appear in various guises, each with subtle
idiosyncrasies warranting particular methodological considerations. This work
considers the intersection of two such sub-types of network data. On the one
hand, we consider when a sample of independent networks is observed. Such
data, often assumed to be drawn from some population of networks, is becoming
increasingly prevalent in the neuroscience literature (Behrens and Sporns,
2012; Chung et al., 2021), and has motivated recent methodological
developments on multiple-network models (Lunagómez et al., 2021; Le et al.,
2018; Durante et al., 2017). On the other hand, we assume that paths or edges
represent the units of observation within each network. Typically referred to
as interaction networks, work has similarly been done towards defining
methodologies suitable to their analysis, most notably with the proposal of
so-called edge-exchangeable models (Crane and Dempsey, 2018; Cai et al., 2016;
Caron and Fox, 2017). As far as we are aware, the intersection of these two
cases, that is, where one observes multiple independent interaction networks,
is yet to be considered in the literature.

As a motivating example, consider the Foursquare check-in data set of Yang et
al. (2015). Foursquare is a location-based social network (LSBN) where users
share places they have visited with their friends by ‘checking-in’ to
locations they visit, such as restaurants or music venues. The data set of
Yang et al. (2015) contains historical check-ins of users in New York and
Tokyo. Considering a single user, we can see a day of check-ins as a path
through the set of venue categories, as illustrated in Figure 1. Over an
extended time period, we expect to observe check-ins on multiple days, leading
to a series of paths being observed. In this way, the data of a single user
can be seen as an interaction network. Furthermore, there are data for not one
but possibly hundreds of users, that is, we have multiple interaction
networks.

With a single observation, we are typically interested in analysing its
structure. When faced with an independent sample of networks the following
more familiar statistical questions are raised

  1. (a)

What is an average network?

  2. (b)

How variable are observations about this average?

  3. (c)

Is there heterogeneity in the observations?

These questions are consistent with those arising more generally in object
orientated data analysis (Marron and Dryden, 2021), which considers the
statistical analysis of populations of complex objects. In this work, we focus
on (a) and (b) in particular.

Figure 1: Data of a single user from the Foursquare data set (Yang et al.,
2015). Going from left to right, the first five subplots show the observed
interactions sequence $\mathcal{S}$, with each path representing a single day
of check-ins, whilst the final plot visualises the aggregate multigraph
$\mathcal{G}_{\mathcal{S}}$, where the thickness of an edge is proportional to
the number of times it appears in $\mathcal{G}_{\mathcal{S}}$.

###  1.1 Related Work

We now review closely related work appearing in the literature. Firstly, there
has been work on models suitable for multiple-network data, where observations
are typically represented via graphs. Here, Lunagómez et al. (2021) construct
models through graph distances, using the Fréchet mean and entropy as notions
of the mean and variance respectively. Le et al. (2018), Peixoto (2018) and
Newman (2018) propose measurement error models which view observed networks as
noisy realisations of an unknown ground truth. Along similar lines, Mantziou
et al. (2021) and Young et al. (2022) have recently extended the measurement
error models to capture heterogeneity, providing model-based approaches to
clustering networks.

Others have considered adapting models originally proposed to analyse a single
network. This includes the latent space model (LSM) (Hoff et al., 2002), which
has been extended by Sweet et al. (2013), who assume a hierarchical model in
which each observation is drawn from an LSM with its own parameter, with these
parameters being linked via a prior, Gollini and Murphy (2016), who assume
observations share the same latent coordinates, and Durante et al. (2017), who
take a non-parametric approach, using a mixture of LSMs combined with
shrinkage priors which induce removal of redundant components and unnecessary
dimensions in latent coordinates. The random dot product graph (RDPG) model
(Young and Scheinerman, 2007) has also been extended. Here Levin et al. (2017)
assume observations are drawn i.i.d. from the same RDPG model, whilst Nielsen
and Witten (2018), Wang et al. (2019) and Arroyo et al. (2021) consider
relaxing this i.i.d. assumption, constructing their models to permit variation
in the RDPG parameters across observations, better capturing heterogeneity.
The exponential random graph model (ERGM) (Holland and Leinhardt, 1981; Frank
and Strauss, 1986) has similarly been extended, where Lehmann and White (2021)
consider a hierarchical model (in similar spirit to Sweet et al., 2013),
whilst Yin et al. (2022) consider a finite mixture of ERGMs. Finally, others
have adapted the stochastic blockmodel (SBM) (Nowicki and Snijders, 2001),
with Sweet et al. (2014) building upon their earlier work (Sweet et al.,
2013), assuming a hierarchical model where each observation is drawn from an
SBM with its own parameterisation, whilst Stanley et al. (2016) and Reyes and
Rodriguez (2016) consider mixtures of SBMs.

There has also been work on hypothesis testing for network-valued data
(Ginestet et al., 2017; Durante and Dunson, 2018; Ghoshdastidar et al., 2020;
Chen et al., 2021), where we note in particular Ginestet et al. (2017) make
use of a distance between graphs, in similar spirit to Lunagómez et al. (2021)
and this present work.

All of these works are connected by a desire to answer standard statistical
questions in the context of network-valued data. However, none consider paths
being the units of observation within each network. Instead, all are designed
to analyse network data represented via graphs. As such, to analyse data which
is truly path-observed, their use would require first aggregating observations
to graphs in a pre-processing step; an operation which will often not be
injective and hence lead to a potential loss of information.

In another direction, there is related work which has considered edge or path-
observed network data. In particular, there has been development of so-called
edge-exchangeable network models (Cai et al., 2016; Crane and Dempsey, 2018;
Williamson, 2016; Ghalebi et al., 2019b, a). Of these we note that only Crane
and Dempsey (2018) allow paths as the observational units. Closely related to
the edge-exchangeable models are those based on exchangeable random measures
(Caron and Fox, 2017; Veitch and Roy, 2015). These two streams of work are
connected in so far as they deviate from more traditional models which are
based upon assumptions of vertex exchangeability, and in doing so produce
graphs which exhibit sparsity and heavy-tailed degree distributions; features
often observed empirically. In another direction, others have considered
models based upon higher-order Markov chains (Scholtes, 2017; Peixoto and
Rosvall, 2017), where in line with this present work Scholtes (2017)
considered paths as the units of observation.

The common theme in all these works is a focus on models which can capture a
particular structure within a single network, such as sparsity, heavy-tailed
degree distributions or high-order dependence in paths or edges. They are not,
however, designed to analyse multiple observations. As such, they can only
provide answers to (a)-(c) via post-hoc analysis of parameters inferred for
each observation.

###  1.2 Summary of Contributions

In the literature, it appears there is a present lack of methods to analyse
samples of interaction networks which fully respect the structure of the data;
either one converts observations to graphs, possibly disregarding information,
or one models each observation individually. We look to address this gap. To
this end, we propose a new Bayesian modelling framework. Namely, through use
of distance metrics between observations, we construct families of models via
location and scale parameters, akin to a Gaussian distribution. The location
parameter, itself an interaction network, admits an interpretation analogous
to the mean, whilst the scale parameter can be seen as a notion of variance or
precision. Conducting inference of these parameters thus provides a reasoned
approach to answering questions (a) and (b). Our methodology is intended to
work with any distance metric of the practitioners choosing, leading to a
flexible framework which can be tailored to suit different questions of
interest. We supplement our methodology with specialised Markov chain Monte
Carlo (MCMC) algorithms which facilitate parameter inference. In particular,
by merging the recently proposed involutive MCMC (iMCMC) framework of
Neklyudov et al. (2020) with the exchange algorithm of Murray et al. (2006),
we are able to not only circumvent issues pertaining to doubly-intractable
posteriors, but also navigate a multi-dimensional discrete parameter space.

The remainder of this paper will be structured as follows. In Section 2, we
provide background details regarding the data structure and notation. Then, in
Section 3 we formally introduce our methodology, stating proposed models and
providing examples of distance metrics. In Section 4, we outline our Bayesian
approach to inference, discussing prior specification, stating our assumed
hierarchical model and detailing our proposed MCMC scheme. In Section 5, we
detail simulation studies undertaken to confirm the efficacy of our
methodology and posterior inference scheme, whilst in Section 6 we illustrate
its applicability via an analysis of the Foursquare check-in data. We finalise
with conclusions and discussion in Section 7.

##  2 Data Representation

Due to the nature of the data, observed paths often arrive in a known order
and, depending on the questions one would like to ask, it may or may not be
desirable to encode this in our representation. For example, consider question
(a) of Section 1. When looking to find an average, do we want to take the
observed order into account, finding an average sequence of paths? Or do we
want to disregard the order information, and instead find an average set of
paths? To cover both situations, we propose two data representations which
will be used within our framework. Namely, what we call interaction sequences
and interaction multisets, covering the ordered and un-ordered cases
respectively.

Adopting the vernacular of network data analysis, we formally denote the
collection of entities under consideration via a vertex set $\mathcal{V}$,
assumed to be some discrete set (typically the set of integers
$\mathcal{V}=\\{1,\dots,V\\}$) where we let $V=|\mathcal{V}|$ denote the
number of vertices. For example, in the Foursquare check-in data (Figure 1),
$\mathcal{V}$ would represent the venue categories. An interaction sequence
will be denoted as follows

$\displaystyle\mathcal{S}=(\mathcal{I}_{1},\dots,\mathcal{I}_{N})$

where the $\mathcal{I}_{i}$ will be referred to as interactions. We consider
the case where interactions are represented via paths, that is

$\displaystyle\mathcal{I}_{i}=(x_{i1},\dots,x_{in_{i}})$

where $x_{ij}\in\mathcal{V}$, and thus $\mathcal{S}$ is a sequence of paths.
Returning to the Foursquare example, $\mathcal{I}_{i}$ would denote a single
day of check-ins, where this user started at a venue of category $x_{i1}$ then
moved to category $x_{i2}$ and so on, with $\mathcal{S}$ denoting all the
observed check-ins of a single user over some fixed time period. Furthermore,
the ordering of interactions reflects the order in which they were observed,
for example, $\mathcal{I}_{1}$ appears before $\mathcal{I}_{2}$ and so on.

When the order of interaction arrival is not of interest, we instead represent
the data via an interaction multiset. A multiset is a set with multiplicities,
that is, a set where elements can occur more than once, and is the natural
order-invariant generalisation of a sequence. An interaction multiset will be
denoted as follows

$\displaystyle\mathcal{E}=\\{\mathcal{I}_{1},\dots,\mathcal{I}_{N}\\},$

where the $\\{\\}$ parenthesis signify this is a multiset, where
$\mathcal{I}_{i}$ similarly denote paths. For example, regarding the
Foursquare data, this would simply represent the collection of observed check-
ins for a single user, with the order of interactions as written above
implying nothing with regards to the order of interacton arival, for example,
$\mathcal{I}_{1}$ was not necessarily observed before $\mathcal{I}_{2}$. We
note this representation has strong similarities with that used in Crane and
Dempsey (2018). In particular, what is defined therein as an interaction
network can be seen as a countably infinite interaction multiset.

In this work, we consider the case where multiple interaction sequences or
multisets are observed. For example, in the Foursquare dataset we have check-
in information on not one but many users, and thus observe a sample of
interaction sequences or multisets, one for each user. Representing the $i$th
observation via a interaction sequence $\mathcal{S}^{(i)}$ or multiset
$\mathcal{E}^{(i)}$, and letting $n$ denote the sample size, we therefore
observe

$\displaystyle\mathcal{S}^{(1)},\dots,\mathcal{S}^{(n)}$ or
$\displaystyle\mathcal{E}^{(1)},\dots,\mathcal{E}^{(n)},$

where the choice of representation depends on ones interest in order. Our
methodology will provide a means to analyse such samples of data.

We finish this subsection by discussing aggregation. Both multisets and
sequences of paths can be aggregated to form graphs, and use of any currently
proposed multiple-network methodology would actually necessitate this as a
pre-processing step. A graph
$\mathcal{G}=(\mathcal{E}_{\mathcal{G}},\mathcal{V}_{\mathcal{G}})$ consists
of a set $\mathcal{V}_{\mathcal{G}}$ of vertices and a set
$\mathcal{E}_{\mathcal{G}}$ of edges where
$e=(i,j)\in\mathcal{E}_{\mathcal{G}}$ if there is an edge from vertex
$i\in\mathcal{V}_{\mathcal{G}}$ to $j\in\mathcal{V}_{\mathcal{G}}$. Graphs can
be un-directed, where
$(i,j)\in\mathcal{E}_{\mathcal{G}}\iff(j,i)\in\mathcal{E}_{\mathcal{G}}$, or
they can be directed, where $(i,j)\in\mathcal{E}_{\mathcal{G}}$ need not imply
$(j,i)\in\mathcal{E}_{\mathcal{G}}$ (and vice versa). Any interaction sequence
$\mathcal{S}$ over vertices $\mathcal{V}$ can be aggregated to form a directed
graph $\mathcal{G}_{\mathcal{S}}=(\mathcal{E}_{\mathcal{S}},\mathcal{V})$ as
follows: let $(i,j)\in\mathcal{E}_{\mathcal{S}}$ if the traversal from vertex
$i$ to vertex $j$ was observed at least once in $\mathcal{S}$, that is, if
$x_{kl}=i$ and $x_{k(l+1)}$ for some $1\leq k\leq N$ and $1\leq l\leq n_{k}$.

Within an interaction sequence, a traversal between two vertices may occur
more than once. Thus we can also aggregate to form a so-called multigraph,
where edges can appear any number of times. We denote a multigraph similarly
via $\mathcal{G}=(\mathcal{E}_{\mathcal{G}},\mathcal{V}_{\mathcal{G}})$, where
in this case $\mathcal{E}_{\mathcal{G}}$ is a multiset of edges, that is, any
given edge $(i,j)\in\mathcal{V}\times\mathcal{V}$ can appear in
$\mathcal{E}_{\mathcal{G}}$ multiple times. Any interaction sequence
$\mathcal{S}$ over vertices $\mathcal{V}$ can be aggregated to form a directed
multigraph $\mathcal{G}_{\mathcal{S}}=(\mathcal{E}_{\mathcal{S}},\mathcal{V})$
by including the edge $(i,j)$ each time the traversal from vertex $i$ to
vertex $j$ is observed, or more formally one can let

$\mathcal{E}_{\mathcal{S}}=\left\\{\left(x_{kl},x_{k(l+1)}\right)\,:\,1\leq
k\leq N,\,1\leq l\leq(n_{k}-1)\right\\},$

for example, an aggregate multigraph for the Foursquare check-in data of a
single user can be seen in Figure 1.

One can also aggregate an interaction mutliset $\mathcal{E}$ over vertex set
$\mathcal{V}$ to form a graph or multigraph in the same manner. Suppose
$\tilde{\mathcal{S}}$ is an interaction sequence obtained by placing the paths
of $\mathcal{E}$ in arbitrary order, then we let

$\mathcal{G}_{\mathcal{E}}=\mathcal{G}_{\tilde{\mathcal{S}}}=(\mathcal{E}_{\tilde{\mathcal{S}}},\mathcal{V}),$

which applies equally to the definition of the graph or multigraph, where in
the former case $\mathcal{E}_{\tilde{\mathcal{S}}}$ will be a set of edges,
whilst in the latter it will be a multiset of edges.

Observe that aggregation is not injective, that is, one may have
$\mathcal{S}\not=\mathcal{S}^{\prime}$ (respectively
$\mathcal{E}\not=\mathcal{E}^{\prime})$ whilst
$\mathcal{G}_{\mathcal{S}}=\mathcal{G}_{\mathcal{S}^{\prime}}$ (respectively
$\mathcal{G}_{\mathcal{E}}=\mathcal{G}_{\mathcal{E}^{\prime}}$). For
interaction sequences, a trivial case would be a re-ordering of paths. This,
however, is not the only example, and there can instead be interaction
sequences or multisets which are structurally dissimilar but nonetheless have
equivalent aggregate graphs. In this way, aggregation incurs a loss of
information that will be avoided with our methodology.

##  3 Metric-Based Interaction-Network Models

We now introduce our proposed models for interaction sequences and multisets
respectively. The core idea behind both is an assumption that observed data
points are ‘noisy’ realisations of some (unknown) ground truth, where
quantification of this noise is facilitated via a pre-specified distance
metric. Equivalently, they can be seen as Gaussian-like distributions over
their respective discrete metric spaces, controlled by a location parameter,
itself an interaction sequence or multiset, and a real-valued scale parameter.

###  3.1 Interaction-Sequence Models

We let $\bm{\mathcal{S}}^{*}$ denote the infinite discrete space consisting of
all interaction sequences over the fixed vertex set $\mathcal{V}$, further
details of which can be found in Appendix A. Towards eliciting a probability
distribution over $\bm{\mathcal{S}}^{*}$, we first endow it with a distance
metric
$d_{S}:\bm{\mathcal{S}}^{*}\times\bm{\mathcal{S}}^{*}\to\mathbb{R}_{+}$,
defining a metric space $(\bm{\mathcal{S}}^{*},d_{S})$. We then select an
element of the space $\mathcal{S}^{m}\in\bm{\mathcal{S}}^{*}$, referred to as
the mode, upon which to center the distribution, before choosing $\gamma>0$,
referred to as the diffusion, controlling the concentration of probability
mass in $\bm{\mathcal{S}}^{*}$ about the mode $\mathcal{S}^{m}$. In this way,
$\mathcal{S}^{m}$ and $\gamma$ can be seen as location and scale parameters,
respectively, analogous to the mean and variance of a Gaussian distribution. A
family of probability distributions can now be defined as follows.

######  Definition 1 (Spherical Interaction Sequence Family)

For a given metric $d_{S}(\cdot,\cdot)$ on $\bm{\mathcal{S}}^{*}$, a mode
$\mathcal{S}^{m}\in\bm{\mathcal{S}}^{*}$ and a diffusion parameter $\gamma>0$,
the probability of observing $\mathcal{S}$ is given by

$p(\mathcal{S}\,|\,\mathcal{S}^{m},\gamma)\propto\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\},$ (1)

and we write

$\mathcal{S}\sim\text{SIS}(\mathcal{S}^{m},\gamma)$

if we assume $\mathcal{S}$ was sampled via (1). This we refer to as the
Spherical Interaction Sequence (SIS) family of probability distributions over
$\bm{\mathcal{S}}^{*}$ with parameters $\mathcal{S}^{m}$ and $\gamma$.

###  3.2 Interaction-Multiset Models

For interaction multisets, the family of models we propose are more-or-less a
carbon copy of the SIS models, albeit with slightly different notation. In
this case, $\bm{\mathcal{E}^{*}}$ will denote the sample space; here
consisting of all interaction multisets over the fixed vertex set
$\mathcal{V}$ (further details provided in Appendix A). Again, we endow
$\bm{\mathcal{E}^{*}}$ with a distance metric
$d_{E}\,:\bm{\mathcal{E}^{*}}\times\bm{\mathcal{E}^{*}}\to\mathbb{R}_{+}$,
similarly defining a metric space $(\bm{\mathcal{E}^{*}},d_{E})$. We then
construct a distribution over $\bm{\mathcal{E}^{*}}$ via location and scale in
exactly the same way, where in this case our location parameter will be an
interaction multiset $\mathcal{E}^{m}\in\bm{\mathcal{E}^{*}}$. The resultant
family of probability distributions is defined as follows.

######  Definition 2 (Spherical Interaction Multiset Family)

Given a metric $d_{E}(\cdot,\cdot)$ on $\bm{\mathcal{E}^{*}}$, a mode
$\mathcal{E}^{m}\in\bm{\mathcal{E}^{*}}$, and a diffusion parameter
$\gamma>0$, the probability of observing $\mathcal{E}$ is given by

$p(\mathcal{E}\,|\,\mathcal{E}^{m},\gamma)\propto\exp\\{-\gamma
d_{E}(\mathcal{E},\mathcal{E}^{m})\\},$ (2)

and we write

$\mathcal{E}\sim\text{SIM}(\mathcal{E}^{m},\gamma)$

if we assume $\mathcal{E}$ was sampled via (2). This we refer to as the
Spherical Interaction Multiset (SIM) family of probability distributions over
$\bm{\mathcal{E}^{*}}$ with parameters $\mathcal{E}^{m}$ and $\gamma$.

Note that in (1) and (2) no reference is made to $d_{S}(\cdot,\cdot)$ or
$d_{E}(\cdot,\cdot)$, even though the respective distributions clearly depend
on them. The reasoning here is these values are not intended to be model
parameters but instead subjective choices made by the practitioner prior to
any analysis.

###  3.3 Properties of Interaction-Network Models

We now highlight some key aspects of both model families. Being almost
identical, these will be shared between the two, and thus for brevity we here
provide details regarding the interaction-sequence models only, noting that
the same properties will hold for the multiset models, up to a change of
notation.

Firstly, distributions of the form (1) agree with the intuition of location
and scale. In particular, for fixed $d_{S}(\cdot,\cdot)$, an observation
$\mathcal{S}$ has the highest probability when $\mathcal{S}=\mathcal{S}^{m}$,
and as $\mathcal{S}$ moves further away from $\mathcal{S}^{m}$, that is, as
$d_{S}(\mathcal{S},\mathcal{S}^{m})$ increases, we see a decay in probability.
Thus, $\mathcal{S}^{m}$ indeed represents the mode of the distribution and can
be seen as controlling the location of probability within
$\bm{\mathcal{S}}^{*}$. Furthermore, the rate at which the probability
decreases with respect to distance is controlled by $\gamma$, with larger
values leading to a faster decrease, implying a greater concentration of
probability mass about the mode. In this way, $\gamma$ is analogous to the
inverse of the variance, often referred to as the precision.

This latter aspect, the control of variance by $\gamma$, can be formalised as
in Lunagómez et al. (2021) via the entropy, defined to be

$\displaystyle H(\mathcal{S}^{m},\gamma)$
$\displaystyle:=-\mathbb{E}\left[\log
p(\mathcal{S}\,|\,\mathcal{S}^{m},\gamma)\right],$ (3)

which can be seen to quantify the uniformity of
$p(\mathcal{S}|\mathcal{S}^{m},\gamma)$, whereby larger values of
$H(\mathcal{S}^{m},\gamma)$ imply this distribution is ‘more uniform’ over
$\bm{\mathcal{S}}^{*}$, with a minimum value of $H(\mathcal{S}^{m},\gamma)=0$
attained by a pointmass. The entropy also has an interpretation with regards
to randomness or variance, whereby distributions with a higher entropy are
more random, that is more variable. It can be shown that with any
$d_{S}(\cdot,\cdot)$ and $\mathcal{S}^{m}$, the entropy
$H(\mathcal{S}^{m},\gamma)$ is a monotonic function of $\gamma$ (Appendix S1).
Consequently, one can say $\gamma$ controls the variance of the distribution,
as formalised via the entropy.

The distribution as stated in Definition 1 is normalised as follows

$p(\mathcal{S}\,|\,\mathcal{S}^{m},\gamma)=Z(\mathcal{S}^{m},\gamma)^{-1}\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}$ (4)

where

$Z(\mathcal{S}^{m},\gamma)=\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}$

is the normalising constant, often referred to as the partition function. In
general, this summation is intractable; an aspect which will come into play
significantly when we consider the computational aspects of the methodology in
later sections. In fact, with $\bm{\mathcal{S}}^{*}$ being infinite, there is
no guarantee that (3.3) will even exist for a given $\gamma$. Consequently,
for some parameterisations (4) may be an improper distribution. Pragmatically,
this can cause issues when sampling from these models, whereby one can observe
a divergence in dimension for low values of $\gamma$. However, this will
generally not be an issue provided one chooses $\gamma$ sufficiently large. In
any case, it can be helpful to constrain the sample space to be finite to
ensure existence of normalising constants, details of which can be found in
Section A.1.

We conclude this section by highlighting similarities of our models with those
seen in other applications. Firstly, as already touched upon, both families
are a direct extension of the spherical network family (SNF) proposed by
Lunagómez et al. (2021) for graphs. Also, our assumption that observed data
points are ‘noisy’ realisations of an unknown ground truth is very much in
keeping with the measurement error models for graphs (Le et al., 2018;
Peixoto, 2018; Newman, 2018). Finally, there are connections with the complex
Watson distribution (Mardia and Dryden, 1999), used in the context of
statistical shape analysis, and the Mallows model (Vitelli et al., 2018) used
for preference learning.

###  3.4 Distance Metrics

A fundamental aspect of our modelling framework is the use of distance
metrics. Not only are these key in facilitating the construction of both model
families but they also represent an important modelling choice. In this
section, examples of two distances measures usable within our framework will
be given for interaction sequences and multisets, respectively.

For interaction sequences, one can consider the so-called edit distance (Bolt
et al., 2022, Definition 4.1). This admits two interpretations: the minimum
cost of transforming one sequence into the other, or the cost of an optimal
pairing of paths across the two sequences. It is the latter interpretation
that we formally define. Suppose $\mathcal{S}$ and $\mathcal{S}^{\prime}$
denote two interaction sequences, then a pairing of their paths can be encoded
via a matching, defined to be a multiset of pairs

$\mathcal{M}=\\{(\mathcal{I},\mathcal{I}^{\prime})\,:\,\mathcal{I}\in\mathcal{S},\mathcal{I}^{\prime}\in\mathcal{S}^{\prime}\\},$

satisfying certain properties. For example, each interaction of one
observation may be paired with at most one of the other and vice versa. For
comparison of sequences, it makes to sense to further constrain the structure
of matchings so as to ensure the preservation of order. This is achieved by
allowing only so-called monotone matchings (see Figure 2).

The edit distance between $\mathcal{S}$ and $\mathcal{S}^{\prime}$ is now
defined as follows

$\displaystyle
d_{\text{Edit}}(\mathcal{S},\mathcal{S}^{\prime}):=\min_{\mathcal{M}\in\bm{\mathcal{M}^{*}}_{m}}\left\\{\sum_{(\mathcal{I},\mathcal{I}^{\prime})\in\mathcal{M}}d_{I}(\mathcal{I},\mathcal{I}^{\prime})+\sum_{\mathcal{I}\in\mathcal{M}^{c}_{\mathcal{S}}}d_{I}(\mathcal{I},\Lambda)+\sum_{\mathcal{I}^{\prime}\in\mathcal{M}^{c}_{\mathcal{S}^{\prime}}}d_{I}(\Lambda,\mathcal{I}^{\prime})\right\\}$
(5)

where $\bm{\mathcal{M}^{*}}_{m}$ is the set of monotone matchings between the
two sequences; $\mathcal{M}_{\mathcal{S}}$ is the restriction of $\mathcal{M}$
to $\mathcal{S}$, that is, the set of interactions in $\mathcal{S}$ included
in the matching, with
$\mathcal{M}^{c}_{\mathcal{S}}:=\mathcal{S}\setminus\mathcal{M}_{\mathcal{S}}$
the interactions of $\mathcal{S}$ not included in $\mathcal{M}$; $d_{I}$ is
some distance metric between interactions, that is, paths; and $\Lambda$
denotes the empty interaction. That is, we sum the pairwise distances of
matched interactions (equivalent to the cost of transforming interactions to
one another), and then sum penalisations for unmatched interactions, where an
unmatched interaction $\mathcal{I}$ incurs penalty
$d_{I}(\mathcal{I},\Lambda)$ (equivalent to the cost of inserting or deleting
$\mathcal{I}$). Computation of $d_{\text{Edit}}$ can be achieved via dynamic
programming (Bolt et al., 2022) at a complexity of $\mathcal{O}(NM)$, where
$N$ and $M$ are the lengths of $\mathcal{S}$ and $\mathcal{S}^{\prime}$
respectively.

For interaction multisets, we consider the matching distance (Bolt et al.,
2022, Definition 3.1). This has close connections with the edit distance, with
the interpretation being more-or-less identical: on the one hand, the minimum
cost of editing one multiset to get the other, whilst on the other an optimal
matching of interactions across the two observations. The differential factor
is an absence of order information, making consideration of monotonicity in
the matching unnecessary. Given two interaction multisets $\mathcal{E}$ and
$\mathcal{E}^{\prime}$, the matching distance is defined as follows

$\displaystyle
d_{\text{M}}(\mathcal{E},\mathcal{E}^{\prime}):=\min_{\mathcal{M}\in\bm{\mathcal{M}^{*}}}\left\\{\sum_{(\mathcal{I},\mathcal{I}^{\prime})\in\mathcal{M}}d_{I}(\mathcal{I},\mathcal{I}^{\prime})++\sum_{\mathcal{I}\in\mathcal{M}^{c}_{\mathcal{E}}}d_{I}(\mathcal{I},\Lambda)+\sum_{\mathcal{I}^{\prime}\in\mathcal{M}^{c}_{\mathcal{E}^{\prime}}}d_{I}(\Lambda,\mathcal{I}^{\prime})\right\\}$
(6)

where $\bm{\mathcal{M}^{*}}$ is the set of matchings between the two multisets
(with no monotonicity constraints), whilst $\mathcal{M}_{\mathcal{E}}$ is the
restriction of $\mathcal{M}$ to $\mathcal{E}$, that is, the subset of
interactions from $\mathcal{E}$ included in $\mathcal{M}$, with
$\mathcal{M}^{c}_{\mathcal{E}}:=\mathcal{E}\setminus\mathcal{M}_{\mathcal{E}}$
again denoting the interactions not included in the matching $\mathcal{M}$.
Computation of $d_{\text{M}}$ can be achieved via the Hungarian algorithm
(Kuhn, 1955), with a complexity of $\mathcal{O}(\max(N,M)^{4}+NM)$, where $N$
and $M$ denote the cardinalities of $\mathcal{E}$ and $\mathcal{E}^{\prime}$
respectively (Bolt et al., 2022).

(a) monotone matching

(b) general matching

Figure 2: A comparison of matchings. In both (a) and (b) we have two
interaction sequences displayed displayed top-down. Here (a) shows a monotone
matching of interactions, where the order is preserved, whilst (b) shows a
general matching, where any two interactions can be paired together.

Use of both $d_{\text{Edit}}$ and $d_{\text{M}}$ requires specification of a
distance $d_{I}$ between paths. Here we provide examples of two such
distances, namely the longest common subpath (LSP) and longest common
subsequence (LCS) distances. The core idea of both is to quantify the
dissimilarity of two paths by considering how much they have (or do not have)
in common.

For $\mathcal{I}=(x_{1},\dots,x_{n})$ a subpath of $\mathcal{I}$ from index
$i$ to $j$ is denoted as follows

$\mathcal{I}_{i:j}=(x_{i},\dots,x_{j})$

where $1\leq i\leq j\leq n$. Given another path
$\mathcal{I}^{\prime}=(y_{1},\dots,y_{m})$ one can consider finding subpaths
that both $\mathcal{I}$ and $\mathcal{I}^{\prime}$ have in common, that is

$\mathcal{I}_{i:j}=\mathcal{I}^{\prime}_{k:l}$

where $1\leq i\leq j\leq n$ and $1\leq k\leq l\leq m$ (Figure 3(a)).
Alternatively, one can consider common subsequences. Assuming that
$\bm{v}=(v_{1},\dots,v_{s})$ with $1\leq v_{1}<v_{2}<\dots<v_{s}\leq n$, then
the subsequence of $\mathcal{I}$ indexed by $\bm{v}$ is given by

$\mathcal{I}_{\bm{v}}=(x_{v_{1}},\dots,x_{v_{s}}),$

and one can similarly consider finding subsequences that both $\mathcal{I}$
and $\mathcal{I}^{\prime}$ have in common, that is

$\mathcal{I}_{\bm{v}}=\mathcal{I}^{\prime}_{\bm{u}}$

where $\bm{u}=(u_{1},\dots,u_{r})$ is such that $1\leq
u_{1}<u_{2}<\dots<u_{r}\leq m$ and $\bm{v}$ is as above (Figure 3(b)).

Towards defining a distance, one can consider finding a maximal subpath or
subsequence, that is, one such that there exists no other subpath or
subsequence of greater length. For example, the common subpaths and
subsequences of Figure 3 are maximal. Intuitively, a distance is now defined
by counting the number of entries not included in the common subpath (or
subsequence) from either $\mathcal{I}$ or $\mathcal{I}^{\prime}$.

More formally, the two distances can be defined as follows

$\displaystyle
d_{\text{LSP}}(\mathcal{I},\mathcal{I}^{\prime}):=n+m-2\delta_{\text{LSP}}$
$\displaystyle
d_{\text{LCS}}(\mathcal{I},\mathcal{I}^{\prime}):=n+m-2\delta_{\text{LCS}}$

where $\delta_{\text{LSP}}$ and $\delta_{\text{LCS}}$ denote the maximum
common subpath and subsequence lengths respectively. This equates to counting
the underlined entries of Figure 3, whereby we consider the $n+m$ total
entries (of both paths), and then remove the double counting made for the
entries of the common subpath or subsequence. As an example, in Figure 3 we
have $\delta_{\text{LSP}}=3$ and $\delta_{\text{LCS}}=4$, and thus
$d_{\text{LSP}}(\mathcal{I},\mathcal{I}^{\prime})=7$ whilst
$d_{\text{LCS}}(\mathcal{I},\mathcal{I}^{\prime})=5$.

Computation of these distances reduces to finding the maximum length common
subpaths and subsequences, both of which can be achieved via dynamic
programming with a complexity $\mathcal{O}(nm)$. Further details on the
implementation can be found in Bolt et al. (2022), Appendix C.

(a) Common subpath

(b) Common subsequence

Figure 3: A comparison of common subsequences and subpaths. In (a) and (b) we
see the same pair of paths, with (a) highlighting a common subpath, as
indicated by shaded (green) entries, whilst (b) shows a common subsequence. In
both cases, these are in fact maximal.

###  3.5 Modelling Choices

A key modelling choice is whether to opt for an SIS or SIM model. This should
be framed around ones interest in order information, whereby if one would like
to take order into account then an SIS model should be assumed, choosing an
SIM model in the converse case. This can be equated with making a form of
exchangeability assumption reminiscent of edge-exchangeable network models
(Cai et al., 2016; Crane and Dempsey, 2018), where in assuming an SIM model
one is considering observations equal up to a permutation of path order as
having the same probability of being sampled.

Another important modelling decision is the choice of distance metric. This
will imply certain assumptions regarding the structure of noise in the model,
in turn altering the interpretations of inference. For example, taking
$d_{I}=d_{\text{LSP}}$ within either the edit or matching distances of Section
3.4 will imply samples drawn from the SIS or SIM models will contain paths
which share common subpaths with those of the mode, whilst if one took
$d_{I}=d_{\text{LCS}}$ then samples would instead share common subsequences
with the mode. As a consequence, when inferring model parameters, in
particular the mode, the interpretation thereof will differ. In particular, in
assuming an LSP distance between paths the inferred mode will likely consist
of a multiset or sequence of subpaths often seen together in the observed
data, whilst if we assume an LCS distance we will instead uncover a multiset
or sequence of often-observed subsequences.

##  4 Bayesian Inference

Given an assumed model, the goal of inference is to discern which parameters
are likely to have generated the observed data. In our case, this amounts to
inferring the mode and dispersion parameters. We approach this task via a
Bayesian perspective, first assuming prior distributions for model parameters
before incorporating observed data to form the posterior. Through a
specialised MCMC algorithm, we subsequently obtain samples from the posterior
upon which to base our inference. In this section, we provide details
regarding each of these aspects.

For brevity, we give details regarding the interaction-sequence models
(Section 3.1) only, noting the approach taken for the multiset models is
almost identical, albeit with a change of notation and some minor alterations
to the MCMC algorithms. Full details of inference for the interaction-multiset
models analogous to those given below can be found in Appendix S4.

###  4.1 Priors, Hierarchical Model and Posterior

In specifying a prior for the mode we follow Lunagómez et al. (2021) and
assume it was itself sampled from an SIS model, that is

$\displaystyle\mathcal{S}^{m}\sim\text{SIS}(\mathcal{S}_{0},\gamma_{0})$ (7)

where $(\mathcal{S}_{0},\gamma_{0})$ are specified hyperparameters. For the
dispersion $\gamma$ we simply require a distribution $p(\gamma)$ whose support
is a subset of the non-negative reals. For example, we typically take
$\gamma\sim\text{Gamma}(\alpha_{0},\beta_{0})$ with $(\alpha_{0},\beta_{0})$
being hyperparameters. Given these specifications, an observed sample
$\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}$ is thus assumed to be drawn via

$\displaystyle\mathcal{S}^{(i)}\,|\,\mathcal{S}^{m},\gamma$
$\displaystyle\sim\text{SIS}(\mathcal{S}^{m},\gamma)$ $\displaystyle(\text{for
}i=1,\dots,n)$ (8) $\displaystyle\mathcal{S}^{m}$
$\displaystyle\sim\text{SIS}(\mathcal{S}_{0},\gamma_{0})$
$\displaystyle\gamma$ $\displaystyle\sim p(\gamma).$

The likelihood of the sample $\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}$ is given by

$\displaystyle
p(\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}\,|\,\mathcal{S}^{m},\gamma)$
$\displaystyle=\prod_{i=1}^{n}p(\mathcal{S}^{(i)}\,|\,\mathcal{S}^{m},\gamma)$
$\displaystyle=Z(\mathcal{S}^{m},\gamma)^{-n}\exp\left\\{-\gamma\sum_{i=1}^{n}d_{S}(\mathcal{S}^{(i)},\mathcal{S}^{m})\right\\},$

and we have the following posterior, up to a constant of proportionality

$\displaystyle
p(\mathcal{S}^{m},\gamma\,|\,\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})$
$\displaystyle\propto
p(\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}\,|\,\mathcal{S}^{m},\gamma)p(\mathcal{S}^{m})p(\gamma)$
(9) $\displaystyle\propto
Z(\mathcal{S}^{m},\gamma)^{-n}\exp\left\\{-\gamma\sum_{i=1}^{n}d_{S}(\mathcal{S}^{(i)},\mathcal{S}^{m})\right\\}$
$\displaystyle\qquad\times\exp\\{-\gamma_{0}d_{S}(\mathcal{S}^{m},\mathcal{S}_{0})\\}p(\gamma).$

###  4.2 Sampling from the Posterior

Figure 4: Summary of our MCMC scheme to sample from the SIS posterior. We
first update the mode via the iExchange algorithm, doing an edit allocation
move with probability $\beta$, or a path insertion and deletion move
otherwise. We then update the dispersion via the exchange algorithm.

To sample from the posterior (9), we use a component-wise MCMC algorithm which
alternates between sampling from the two conditional distributions

$\displaystyle
p(\mathcal{S}^{m}\,|\,\gamma,\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})$ and
$\displaystyle
p(\gamma\,|\,\mathcal{S}^{m},\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}).$ (10)

Since the normalising constant $Z(\mathcal{S}^{m},\gamma)$ depends on the
parameters of interest this implies (9) is doubly intractable (Murray et al.,
2006). Such terms will also persist in both conditionals above, making them
also doubly intractable. This precludes the use of standard MCMC algorithms
such as Metropolis-Hastings and necessitates use of the exchange algorithm
proposed by Murray et al. (2006).

A high-level summary of our scheme is visualised in Figure 4. For the
dispersion conditional, being a distribution over the real line, we can apply
the exchange algorithm directly. In contrast, the mode $\mathcal{S}^{m}$ is a
discrete object, the dimensions of which can vary both in terms of the number
of paths and their lengths. This makes the sample space for the mode
conditional far less trivial, and so we consider merging the exchange
algorithm with the involutive MCMC (iMCMC) framework of Neklyudov et al.
(2020); defining what we call the iExchange algorithm. To fully explore the
sample space, we mix together two iExchange moves. In particular, with
probability $\beta$, we enact a move perturbing the paths currently present,
whilst with probability $(1-\beta)$ we attempt a move which varies the number
of paths.

###  4.3 Updating the Dispersion

Here we outline our MCMC scheme to sample from the dispersion conditional. In
this instance, we suppose $\mathcal{S}^{m}$ is fixed and
$q(\gamma^{\prime}|\gamma)$ is some proposal density. In a single iteration,
given current state $\gamma$ we do the following

  1. 1.

Sample a proposal $\gamma^{\prime}$ from $q(\gamma^{\prime}|\gamma)$

  2. 2.

Sample an auxiliary dataset $\\{\mathcal{S}^{*}_{i}\\}_{i=1}^{n}$ of size $n$
(same as observed data) where

$\mathcal{S}_{i}^{*}\overset{\text{i.i.d.}}{\sim}\text{SIS}(\mathcal{S}^{m},\gamma^{\prime}),$

  3. 3.

Evaluate the following probability

$\displaystyle\alpha(\gamma,\gamma^{\prime})$
$\displaystyle=\min\left\\{1,\frac{p(\gamma^{\prime}|\mathcal{S}^{m},\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})p(\\{\mathcal{S}^{*}_{i}\\}_{i=1}^{n}|\mathcal{S}^{m},\gamma)q(\gamma|\gamma^{\prime})}{p(\gamma|\mathcal{S}^{m},\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})p(\\{\mathcal{S}^{*}_{i}\\}_{i=1}^{n}|\mathcal{S}^{m},\gamma^{\prime})q(\gamma^{\prime}|\gamma)}\right\\}$
(11)

  4. 4.

Move to state $\gamma^{\prime}$ with probability
$\alpha(\gamma,\gamma^{\prime})$, staying at $\gamma$ otherwise.

For the proposal $q(\gamma^{\prime}|\gamma)$ we consider sampling
$\gamma^{\prime}$ uniformly over a $\varepsilon$-neighbourhood of $\gamma$
with reflection at zero. More specifically, we first sample
$\gamma^{*}\sim\text{Uniform}(\gamma-\varepsilon,\gamma+\varepsilon)$ and then
let $\gamma^{\prime}=\gamma^{*}$ if $\gamma^{*}>0$ and let
$\gamma^{\prime}=-\gamma^{*}$ otherwise.

This is a direct application of the exchange algorithm (Murray et al., 2006)
and as such the resultant Markov chain admits
$p(\gamma\,|\,\mathcal{S}^{m},\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})$ as its
stationary distribution. Moreover, this is what one might call an “exact-
approximate” MCMC algorithm, in the sense that (asymptotically) samples drawn
thereof will be distributed according to the desired target, meaning that one
could in theory obtain exact samples given infinite resource. A closed form of
(11) and derivation thereof can be found in Section S3.1.

###  4.4 Updating the Mode

We now outline our MCMC scheme to sample from the mode conditional. The key
difference here is in the proposal generation mechanism, which follows the
iMCMC algorithm (Neklyudov et al., 2020) in using a combination of random
sampling and deterministic maps. Here we assume the dispersion $\gamma$ is
fixed and $\mathcal{S}^{m}$ denotes our current state. Instead of specifying a
proposal density, one defines auxiliary variables $u\in\mathcal{U}$, a
deterministic function
$f\,:\,\bm{\mathcal{S}}^{*}\times\mathcal{U}\to\bm{\mathcal{S}}^{*}\times\mathcal{U}$
and a conditional distribution $q(u\,|\,\mathcal{S}^{m})$ over auxiliary
variables. The function $f$ must also be an involution, meaning that is acts
as its own inverse, that is, $f^{-1}=f$. A single iteration now consists of
the following

  1. 1.

Sample $u\sim q(u|\mathcal{S}^{m})$

  2. 2.

Invoke involution
$f(\mathcal{S}^{m},u)=([\mathcal{S}^{m}]^{\prime},u^{\prime})$, obtaining
proposal $[\mathcal{S}^{m}]^{\prime}$

  3. 3.

Sample auxiliary dataset $\\{\mathcal{S}^{*}_{i}\\}_{i=1}^{n}$ of size $n$
where

$\mathcal{S}^{*}_{i}\overset{\text{i.i.d.}}{\sim}\text{SIS}([\mathcal{S}^{m}]^{\prime},\gamma)$

  4. 4.

Evaluate the following probability

$\displaystyle\alpha(\mathcal{S}^{m},[\mathcal{S}^{m}]^{\prime})$
$\displaystyle=\min\left\\{1,\frac{p([\mathcal{S}^{m}]^{\prime}\,|\,\gamma,\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})p(\\{\mathcal{S}^{*}_{i}\\}_{i=1}^{n}\,|\,\mathcal{S}^{m},\gamma)q(u^{\prime}\,|\,[\mathcal{S}^{m}]^{\prime})}{p(\mathcal{S}^{m}\,|\,\gamma,\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})p(\\{\mathcal{S}^{*}_{i}\\}_{i=1}^{n}\,|\,[\mathcal{S}^{m}]^{\prime},\gamma)q(u\,|\,\mathcal{S}^{m})}\right\\}$
(12)

  5. 5.

Move to state $[\mathcal{S}^{m}]^{\prime}$ with probability
$\alpha(\mathcal{S}^{m},[\mathcal{S}^{m}]^{\prime})$, staying at
$\mathcal{S}^{m}$ otherwise.

Much like the proposal density of a Metropolis-Hasting or exchange algorithm,
the $u$, $f(\mathcal{S}^{m},u)$ and $q(u\,|\,\mathcal{S}^{m})$ represent free
choices. We consider mixing together two such specifications, details of which
we provide in the next section.

This scheme represents an instance of what we call the iExchange algorithm
(Algorithm 1, Appendix S2). As shown in Appendix S2, this can be seen as a
special case of the iMCMC algorithm. As such, this represents an exact-
approximate MCMC algorithm with the resultant Markov chain admitting
$p(\mathcal{S}^{m}|\gamma,\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})$ as its
stationary distribution. Note the iExchange algorithm as defined in Appendix
S2 includes a Jacobian term in the acceptance probability which we do not
include above. The reasoning being that since both $\bm{\mathcal{S}}^{*}$ and
$\mathcal{U}$ are discrete spaces and $f(\mathcal{S},u)$ is a one-to-one
function (since it is invertible) such terms are not required.

###  4.5 Mode Update Moves

We now give details regarding two iExchange specifications for the mode
conditional updates. In the first, we keep the number of paths fixed, varying
only the path lengths or what we call the inner dimension. For example, in the
context of the Foursquare data, this would amount to altering a particular
sequence of check-ins. In the second, we look to vary the number of paths or
what we call the outer dimension. For example, in the Foursquare data this
would equate to introducing or removing a whole day of check-ins.

####  4.5.1 Edit Allocation

Figure 5: Illustrating the edit allocation move. Shaded entries indicate
deletions and insertions, whilst bars visualise allocation of edits to paths.
Bar height is proportional to the number of edits allocated to a path $z$,
whilst the green (top) portion of the bar denotes the number of insertions $a$
and the pink (bottom) portion represents the number of deletions $d$.

Supposing $\mathcal{S}^{m}=(\mathcal{I}_{1},\dots,\mathcal{I}_{N})$ is our
current state, the main idea of this move is to allocate a number of “edits”
to each path in $\mathcal{S}^{m}$. These edits consist of inserting and
deleting entries, where if the number of insertions and deletions are
unbalanced, paths of smaller or larger sizes relative to the current state
will be proposed, thus varying the inner dimension. For an illustration, see
Figure 5.

We now give descriptive details of this proposal generation mechanism and show
how it can be cast in the light of iMCMC. First, we specify the total number
of edits to be made, denoting this $\delta\in\mathbb{Z}_{\geq 1}$. Next, we
specify an allocation of these edits to the paths of $\mathcal{S}^{m}$,
denoting this $\bm{z}=(z_{1},\dots,z_{N})$, where $z_{i}\in\mathbb{Z}_{\geq
0}$ denotes the number of edits allocated to the $i$th path such that
$\sum_{i=1}^{N}z_{i}=\delta$. For example, in Figure 5 we have $\delta=7$ and
$\bm{z}=(3,1,3,0)$.

Given $z_{i}$ we edit the $i$th path $\mathcal{I}_{i}$ to obtain a
corresponding proposal $\mathcal{I}_{i}^{\prime}$ in the following manner.
First, we partition the $z_{i}$ edits between deletions and insertions,
letting $d_{i}\in\\{0,\dots,\min(n_{i},z_{i})\\}$ denote the number of
deletions, where $n_{i}$ denotes the length of the $i$th path, with
$a_{i}=z_{i}-d_{i}$ then denoting the number of insertions. Note, we cannot
delete more entries than are present, hence the restriction
$d_{i}\leq\min(n_{i},z_{i})$.

The penultimate step is to specify which entries to delete and where to insert
new entries, which we denote via subsequences. Introducing the notation
$[n]=(1,\dots,n)$, we define subsequence of $[n]$ of size $m$ to be a vector
$\bm{v}=(v_{1},\dots,v_{m})$ such that $1\leq v_{1}<v_{2}<\dots<v_{m}\leq n$.
Now, we let $\bm{v}_{i}$ be a subsequence of $[n_{i}]$ of size $d_{i}$
denoting the entries of $\mathcal{I}_{i}$ to be deleted, whilst
$\bm{v}^{\prime}_{i}$ is subsequence of $[m_{i}]$ of size $a_{i}$, denoting
the location of entries to be inserted in $\mathcal{I}^{\prime}_{i}$, where
$m_{i}=n_{i}-d_{i}+a_{i}$ denotes the length of $\mathcal{I}^{\prime}_{i}$.
For example, considering the first path in Figure 5 we have
$\mathcal{I}_{1}=(a,b,a,b,a)$ and $\mathcal{I}^{\prime}_{1}=(d,a,b,a)$ with
$\bm{v}_{1}=(4,5)$ and $\bm{v}^{\prime}_{1}=(1)$ indexing the deletions and
insertions respectively. The final step is to specify entries to insert, which
we denote $\bm{y}_{i}=(y_{i1},\dots,y_{ia_{i}})$ where $y_{ij}\in\mathcal{V}$.
For example, in Figure 5 we have $\bm{y}_{1}=(d)$.

Given the information above, one can enact the specified deletions and
insertions, mapping to a proposal
$[\mathcal{S}^{m}]^{\prime}=(\mathcal{I}^{\prime}_{1},\dots,\mathcal{I}^{\prime}_{N})$.
This can be viewed in the iMCMC framework as follows. First, collate all this
information into the auxiliary variable $u=(\delta,\bm{z},u_{1},\dots,u_{N})$
where $u_{i}=(d_{i},\bm{v}_{i},\bm{v}^{\prime}_{i},\bm{y}_{i})$. Now, if we
write the required involution as follows

$f(\mathcal{S}^{m},u)=(f_{1}(\mathcal{S}^{m},u),f_{2}(\mathcal{S}^{m},u))=([\mathcal{S}^{m}]^{\prime},u^{\prime}),$

then in enacting the specified edit operations we have effectively defined the
first component $f_{1}(\mathcal{S}^{m},u)=[\mathcal{S}^{m}]^{\prime}$.
Specification of the second component is more involved, and so we delegate
these details to Section S3.3. Regarding the auxiliary distribution
$q(u|\mathcal{S}^{m})$, we consider the following

$\displaystyle\delta$
$\displaystyle\sim\text{Uniform}\\{1,\dots,\nu_{\text{ed}}\\}$
$\displaystyle\bm{z}\,|\,\delta$
$\displaystyle\sim\text{Multinomial}(\delta\,;\,1/N,\dots,1/N)$ $\displaystyle
d_{i}\,|\,z_{i}$
$\displaystyle\sim\text{Uniform}\\{0,\dots,\min(z_{i},n_{i})\\}$
$\displaystyle(\text{for }i=1,\dots,N)$

whilst $\bm{v}_{i}$ and $\bm{v}^{\prime}_{i}$ are drawn uniformly and the
entry insertions $\bm{y}_{i}$ are assumed to be sampled from some general
distribution $q(\bm{y}_{i}|\mathcal{I}_{i})$, which we typically take to be
the uniform distribution over $\mathcal{V}$. The only tuning parameter here is
$\nu_{\text{ed}}$, which controls the aggressiveness of proposals, with larger
values leading to more edits being attempted on average.

Further details, including full definition of the involution $f$, examples of
possible insertion distributions $q(\bm{y}_{i}|\mathcal{I}_{i})$ and
derivations of key terms appearing the acceptance probability (12), can be
found in Section S3.3.

####  4.5.2 Path Insertion and Deletion

Figure 6: Illustrating path insertion and deletion move, where given current
state $\mathcal{S}^{m}$ the proposal $[\mathcal{S}^{m}]^{\prime}$ is obtained
by deleting and inserting the highlighted paths.

With this move we look to vary the outer dimension, that is, the number of
paths. Similar to Section 4.5.1, we consider doing so by random deletion and
insertion. The difference in this case is that we delete and insert whole
paths (see Figure 6).

In particular, with $\mathcal{S}^{m}=(\mathcal{I}_{1},\dots,\mathcal{I}_{N})$
denoting our current state, we first choose a total number of insertions and
deletions $\varepsilon\in\mathbb{Z}_{\geq 1}$. Next, we partition these,
letting $d\in\\{0,\dots,\min(N,\varepsilon)\\}$ denote the number of
deletions, leaving $a=\varepsilon-d$ insertions. For example, in Figure 6 we
have $\varepsilon=3$, $d=2$ and $a=1$. Next, we choose locations of deletions
and insertions. In particular, we let $\bm{v}$ be a length $d$ subsequence of
$[N]$ denoting which paths of $\mathcal{S}^{m}$ are to be deleted, whilst
$\bm{v}^{\prime}$ is a length $a$ subsequence of $[M]$, where $M=N-d+a$,
denoting where inserted paths will be located in our proposal
$[\mathcal{S}^{m}]^{\prime}$. For example, in Figure 6 we have $\bm{v}=(2,4)$
and $\bm{v}^{\prime}=(3)$. Finally, for each $i=1,\dots,a$ we choose some path
$\mathcal{I}^{*}_{i}$ to insert into entry $v^{\prime}_{i}$ of
$[\mathcal{S}^{m}]^{\prime}$. For example, in Figure 6 we have a single path
$\mathcal{I}^{*}_{1}=(c,b,b,a)$ which we insert to the third entry.

As in Section 4.5.1, given the information above we can insert and delete the
corresponding paths to obtain a proposal $[\mathcal{S}^{m}]^{\prime}$.
Collating this into the auxiliary variable
$u=(\varepsilon,d,\bm{v},\bm{v}^{\prime},\mathcal{I}^{*}_{1},\dots,\mathcal{I}^{*}_{a})$
this can similarly be seen as defining the first component of the required
involution, with details of the second component found in Section S3.4.
Regarding sampling of auxiliary variables we consider the following

$\displaystyle\varepsilon$
$\displaystyle\sim\text{Uniform}\\{1,\dots,\nu_{\text{td}}\\}$ $\displaystyle
d\,|\,\varepsilon$
$\displaystyle\sim\text{Uniform}\\{0,\dots,\min(N,\varepsilon)\\}$

whilst we sample $\bm{v}$ and $\bm{v}^{\prime}$ uniformly and assume path
insertions $\mathcal{I}^{*}_{i}$ are drawn from some general distribution over
paths $q(\mathcal{I}|\mathcal{S}^{m})$. This leaves two tuning parameters,
$\nu_{\text{td}}$ and $q(\mathcal{I}|\mathcal{S}^{m})$, which in combination
facilitate control over the aggressiveness of proposals. In particular,
$\nu_{\text{td}}$ controls the number of deletions and insertions attempted,
whilst $q(\mathcal{I}|\mathcal{S}^{m})$ affects how impactful each of these
insertions and deletions are. Again, further details can be found in Section
S3.4.

###  4.6 Sampling Auxiliary Data

Both algorithms to target the conditionals outlined in Sections 4.3 and 4.4
require exact sampling of auxiliary data from appropriate interaction-sequence
models. Unfortunately, we cannot do this in general. Instead, we consider
replacing this with approximate samples obtained via an iMCMC algorithm.

In particular, suppose we would like to obtain samples from an
$\text{SIS}(\mathcal{S}^{m},\gamma)$ model. Assuming that $\mathcal{S}$
denotes the current state, and auxiliary variables $u$, involution
$f(\mathcal{S},u)$ and auxiliary distribution $q(u|\mathcal{S})$ have be
defined, in a single iteration we do the following

  1. 1.

Sample $u\sim q(u|\mathcal{S})$

  2. 2.

Invoke involution $f(\mathcal{S},u)=(\mathcal{S}^{\prime},u^{\prime})$

  3. 3.

Evaluate the following probability

$\displaystyle\alpha(\mathcal{S},\mathcal{S}^{\prime})$
$\displaystyle=\min\left\\{1,\frac{p(\mathcal{S}^{\prime}\,|\,\mathcal{S}^{m},\gamma)q(u^{\prime}\,|\,\mathcal{S}^{\prime})}{p(\mathcal{S}^{\,}|\,\mathcal{S}^{m},\gamma)q(u\,|\,\mathcal{S})}\right\\}$
(13)

  4. 4.

Move to state $\mathcal{S}^{\prime}$ with probability
$\alpha(\mathcal{S},\mathcal{S}^{\prime})$, staying at $\mathcal{S}$
otherwise.

where $p(\mathcal{S}|\mathcal{S}^{m},\gamma)$ denotes the likelihood as given
in (1). Towards specifying $u$, $f(u,\mathcal{S})$ and $q(u|\mathcal{S})$, we
now recycle the moves of Section 4.4, again mixing these together with some
proportion $\beta\in(0,1)$. Note, as in Section 4.4, we omit the Jacobian term
in the acceptance probability above since we are working with discrete spaces.

In sampling auxiliary data in this manner, we now have two MCMC-based
elements: what one might call the outer MCMC algorithm, navigating the
parameter space, and the inner MCMC algorithm, sampling auxiliary data. We
note this approach has been considered by others. In particular, Liang (2010)
proposed the so-called double Metropolis-Hastings algorithm which replaces the
exact samples of the exchange algorithm with those obtained via a Metropolis-
Hasting scheme. The difference in our case is use of the more general iMCMC
framework, be that in the outer MCMC scheme (as in the iExchange algorithm),
or the inner MCMC scheme (as outlined above).

A consequence of using approximate auxiliary samples within the algorithms of
Sections 4.3 and 4.4 is the resulting schemes will become approximate, as
opposed to exact-approximate. That is to say, even in the theoretical limit,
samples will not necessarily be distributed according to the desired target
but instead an approximation thereof. However, as the auxiliary samples look
more like an i.i.d. sample one will get closer to the respective exact-
approximate algorithm. Thus, one can in theory get arbitrarily close to an
exact-approximate scheme by taking steps to reduce the bias of the MCMC-based
auxiliary samples, such as introducing a burn-in period or taking a lag
between samples.

##  5 Simulation Studies

In this section we outline simulation studies undertaken to confirm the
efficacy of our methodology. In the first two, we examine the posterior
concentration, exploring how this is affected by variability of observed data
and structural features of the mode. In the third, we assess convergence of
the posterior predictive via a missing data problem. In each, we will be
working with the interaction-sequence models.

###  5.1 Posterior Concentration

Figure 7: A summary of our first simulation (Section 5.1), where the top plot
visualises the scale of the SIS model used therein, in particular, for
different values of $\gamma$ it shows
$\\{d_{S}(\mathcal{S}^{(i)},\mathcal{S}^{m}_{\text{true}})\\}_{i=1}^{1000}$
where $\mathcal{S}^{(i)}\sim\text{SIS}(\mathcal{S}^{m}_{\text{true}},\gamma)$,
sampled via the iMCMC scheme of Section 4.6. The remaining two plots summarise
simulation outputs for each pair $(\gamma_{\text{true}},n)$, where the middle
shows distributions of $\bar{d}$, the average distance to the true mode,
whilst the bottom shows $(\bar{\gamma}-\gamma_{\text{true}})$, the bias of the
dispersion posterior mean relative to the truth.

If the data were generated by an SIS model at known parameters, then ideally
the posterior (9) should concentrate about these as the sample size grows. The
goal of the first two simulations is to confirm empirically that such
behaviour is observed and to explore what factors might come into play. The
high-level approach is the following. Given true mode
$\mathcal{S}^{m}_{\text{true}}$ and dispersion $\gamma_{\text{true}}$, we draw
a sample $\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}$ where

$\mathcal{S}^{(i)}\sim\text{SIS}(\mathcal{S}^{m}_{\text{true}},\gamma_{\text{true}})$

before obtaining samples $\\{(\mathcal{S}_{i}^{m},\gamma_{i})\\}_{i=1}^{m}$
from the posterior
$p(\mathcal{S}^{m},\gamma|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})$. We then assess
the behaviour of these samples via the following summary measures

$\displaystyle\bar{d}:=\frac{1}{m}\sum_{i=1}^{m}d_{S}(\mathcal{S}^{m}_{i},\mathcal{S}^{m}_{\text{true}})$
$\displaystyle\bar{\gamma}:=\frac{1}{m}\sum_{i=1}^{m}\gamma_{i}$

where ideally $\bar{d}$ should be close to zero and
$\bar{\gamma}\approx\gamma_{\text{true}}$. By repeating this a number of times
for different $n$ and evaluating these summaries we can thus get a sense of
how the posterior is concentrating about the true parameters.

Now, recall the dispersion works inversely to the variance, in that lower
values lead to more variable data (Figure 7, top). Intuitively, when the data
is more variable it will be harder to discern the true mode
$\mathcal{S}^{m}_{\text{true}}$, and thus we expect $\bar{d}$ to decrease more
slowly for lower values of $\gamma_{\text{true}}$. Alternatively, as can be
seen in Figure 7, when $\gamma_{\text{true}}$ is smaller the difference of
their parameterised distributions (as described by the distribution of
distances to the mode) becomes more marked relative to neighbouring values. As
such, we might also expect smaller values for the dispersion to be easier to
recover.

To explore for such properties, we varied $\gamma_{\text{true}}$ and $n$
whilst keeping $\mathcal{S}^{m}_{\text{true}}$ fixed. In particular, we
considered $\gamma_{\text{true}}=3.5,3.7,4.0,4.3,4.6,4.9$ (highlighted in
Figure 7, top) and $n=25,50,75,100$. The distance we took to be
$d_{S}=d_{\text{Edit}}$ with $d_{I}=d_{\text{LCS}}$ between paths. We fixed
$V=20$, and constrained the sample space as discussed in Section A.1, assuming
at most $L=20$ paths in any observation, with each path being at most length
$K=10$.

The mode $\mathcal{S}^{m}_{\text{true}}$ of length $N=10$ we fixed throughout,
sampled from the Hollywood model of Crane and Dempsey (2018). In particular,
we drew $\mathcal{S}^{m}_{\text{true}}\sim\text{Hollywood}(\alpha,\theta,\nu)$
where

$\displaystyle\alpha=-0.3$ $\displaystyle\theta=0.3V$
$\displaystyle\nu=\text{TrPoisson}(3,1,K),$

where $\text{TrPoisson}(\lambda,a,b)$ denotes a truncated Poisson distribution
with $\lambda>0$ the parameter of a standard Poisson, whilst $0\leq
a<b\leq\infty$ are the lower and upper bounds. This set-up for the Hollywood
model, with $\alpha<0$ and $\theta=-V\alpha$, corresponds to the finite
setting, implying the sampled interaction sequences will have at most $V$
vertices.

Regarding priors, we considered an uninformative set-up with
$(\mathcal{S}_{0},\gamma_{0})=(\hat{\mathcal{S}},0.1)$ where

$\hat{\mathcal{S}}:=\operatorname*{argmin}_{\mathcal{S}\in\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}}\sum_{i=1}^{n}d^{2}_{S}(\mathcal{S}^{(i)},\mathcal{S})$

denotes the sample Fréchet mean, whilst we took
$\gamma\sim\text{Uniform}(0.5,7.0)$. Here we note the sample
$\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}$ used to obtain $\hat{\mathcal{S}}$ will be
different in each repetition of the simulation, and consequently so will
$\hat{\mathcal{S}}$.

Now, for each pair $(\gamma_{\mathrm{true}},n)$ we (i) sampled $n$
observations from an
$\text{SIS}(\mathcal{S}^{m}_{\mathrm{true}},\gamma_{\mathrm{true}})$ model,
using the iMCMC scheme outlined in Section 4.6, with a burn-in period of
50,000 and taking a lag of 500 between samples (ii) obtained $m=250$ samples
from the posterior using the component-wise MCMC scheme of Section 4.2, with a
burn-in period of 25,000 and taking a lag of 100 between samples111One must
also parameterise the MCMC algorithm used to sample the auxiliary data. These
were tuned by considering acceptance probabilities observed when sampling from
an $\text{SIS}(\mathcal{S}^{m}_{\mathrm{true}},\gamma_{\mathrm{true}})$
distribution. (iii) evaluated summary measures $\overline{d}$ and
$\bar{\gamma}$.

We repeated (i)-(iii) 100 times in each case, the results of which are
summarised in Figure 7. Consulting the middle plot, we observe that $\bar{d}$
decreases with $n$ across all cases, indicating a concentration of the
posterior about the true mode. Furthermore, this decrease is more gradual for
lower values of $\gamma_{\text{true}}$, agreeing with intuition. Turning to
the bottom plot, the most obvious feature is bias in $\bar{\gamma}$ relative
to the truth. Note this is expected, since we have used approximate MCMC
samples within our component-wise scheme of Section 4.2. We do, however, see a
reduction in this bias as the sample size grows. Furthermore, for the larger
values of $n$ we begin to see a clearer difference in the variance of
$\bar{\gamma}$ across different values of $\gamma_{\text{true}}$. In
particular, the variance appears to be smaller for lower values of
$\gamma_{\text{true}}$, agreeing with the intuition that these are easier to
estimate.

###  5.2 Effect of Mode Structure

Here we explored whether structural features of the mode might impact its
inference. Adopting the same modelling set-up as the previous simulation, but
in this case fixing the true dispersion to $\gamma_{\text{true}}=4.5$, we re-
sampled the mode in each repetition via

$\mathcal{S}_{\text{true}}\sim\text{Hollywood}(\alpha,-\alpha V,\nu)$

where we again take $V=20$ and $\nu=\text{TrPoisson}(3,1,K)$, whilst
$\alpha<0$.

The key idea underlying the Hollywood model is a ‘rich get richer’ assumption
made when sampling vertices. This results in $\alpha$ admitting an
interpretation regarding the heavy-tailed nature of vertex counts. In
particular, for a given interaction sequence $\mathcal{S}$ and vertex
$v\in\mathcal{V}$ one can define an analogue of the vertex degree (often
defined for graphs) as follows

$k_{\mathcal{S}}(v):=\text{\\# times $v$ appears in $\mathcal{S}$},$

which thus implies, for each $\mathcal{S}$, a sample
$\\{k_{\mathcal{S}}(v)\,:v\in\mathcal{V},\,k_{\mathcal{S}}(v)>0\\}$, similar
in spirit to the degree distribution. Now, $\alpha$ can be seen to control the
heavy-tailedness of this distribution (see Figure 9), whereby when $\alpha$ is
low one tends to see vertices appearing a similar number of times, whilst when
$\alpha$ is larger these counts become disproportionately focused on a smaller
subset of vertices.

In this simulation, we considered $\alpha=-\tilde{\alpha}$ where
$\tilde{\alpha}=1.35,0.75,0.35,0.12,0.06,0.03,0.01$ (details on how these were
chosen can be found in Section B.1) and $n=25,50,75,100$. For each pair
$(\alpha,n)$ in a single repetition we (i) sampled
$\mathcal{S}_{\text{true}}\sim\text{Hollywood}(\alpha,-\alpha V,\nu)$, (ii)
sampled $n$ observations from an
$\text{SIS}(\mathcal{S}_{\text{true}},\gamma_{\text{true}})$ model (iii)
obtained $m=250$ samples from the posterior, and (iv) evaluated summaries. For
(ii) and (iii) we used exactly the same MCMC set-up as in the previous
simulation.

Figure 9 summarises the output of 100 repetitions for each pair $(\alpha,n)$.
For each $\alpha$, we see values for $\bar{d}$ closer to zero as $n$ grows,
indicating concentration about the true mode. Furthermore, $\alpha$ shows no
clear sign of impacting this concentration. Regarding the dispersion posterior
mean $\bar{\gamma}$, as in the previous simulation we observe bias relative to
the truth, with this bias reducing as $n$ grows. Furthermore, this is the same
across all $\alpha$, with no clear sign that $\alpha$ affects the inference of
these values.

Figure 8: Visualising the role of $\alpha$ in the Hollywood model. Each plot
shows an aggregate multigraph $\mathcal{G}_{\mathcal{S}}$ where
$\mathcal{S}\sim\text{Hollywood}(\alpha,-\alpha V,\nu)$ with $V=10$,
$\nu=\text{TrPoisson}(3,1,10)$ and $\alpha$ varying. Edge thickness reflects
edge multiplicity, whilst vertex size is proportional to $k_{\mathcal{S}}(v)$.

Figure 9: Summary of our second simulation (Section 5.2), where for each pair
$(\alpha,n)$ the top subplot shows the distribution of $\bar{d}$, the average
distance to the true mode, whilst the bottom shows the distribution of
$\bar{\gamma}$, the posterior mean dispersion.

###  5.3 Posterior Predictive Efficacy

A desirable feature of the posterior predictive is a growing resemblance of
the true data generating distribution as the sample size increases. In this
simulation, we considered exploring such behaviour in the context of a missing
data problem.

Suppose we have an observation $\mathcal{S}$ in which a single entry is
missing, for example

$\mathcal{S}=((1,2,1,\bullet),(2,3,4,3),(1,2,2,1,2,3))$

with $\bullet$ denoting the unknown entry. Towards predicting its value, let
$\mathcal{S}_{x}$ denote the observation obtained by taking this entry to be
$x\in\mathcal{V}$, that is

$\mathcal{S}_{x}=((1,2,1,x),(2,3,4,3),(1,2,2,1,2,3)),$

and consider assigning a probability to each $x\in\mathcal{V}$ of being the
true entry. If one knew $\mathcal{S}\sim\text{SIS}(\mathcal{S}^{m},\gamma)$,
then such a distribution could be obtained by comparing the relative
probability of $\mathcal{S}_{x}$ for each $x\in\mathcal{V}$, in particular we
could consider

$\displaystyle
p(x|\mathcal{S}^{m},\gamma,\mathcal{S}_{-x}):=\frac{1}{Z(\mathcal{S}^{m},\gamma,\mathcal{S}_{-x})}\exp\\{-\gamma\phi(d_{S}(\mathcal{S}_{x},\mathcal{S}^{m}))\\}$

with
$Z(\mathcal{S}^{m},\gamma,\mathcal{S}_{-x})=\sum_{x\in\mathcal{V}}\exp\\{-\gamma
d_{S}(\mathcal{S}_{x},\mathcal{S}^{m})\\}$ the normalising constant, where we
introduce the notation $\mathcal{S}_{-x}$ to indicate that we are conditioning
on the other known entries (and implicitly also on the dimensions of the
observation). We refer to this as the true predictive for $x\in\mathcal{V}$.

In practice, with the true distribution unknown, one can instead leverage an
observed sample $\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}$ by averaging with respect
to the posterior as follows

$\displaystyle p(x|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n},\mathcal{S}_{-x})$
$\displaystyle=\sum_{\mathcal{S}^{m}\in\bm{\mathcal{S}}^{*}}\int_{\mathbb{R}_{+}}p(x|\mathcal{S}^{m},\gamma,\mathcal{S}_{-x})p(\mathcal{S}^{m},\gamma|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})d\gamma,$

defining the posterior predictive for $x\in\mathcal{V}$, which itself can be
approximated using a sample $\\{(\mathcal{S}^{m}_{i},\gamma_{i})\\}_{i=1}^{m}$
from the posterior via

$\displaystyle\hat{p}(x|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n},\mathcal{S}_{-x})$
$\displaystyle:=\frac{1}{m}\sum_{i=1}^{m}p(x|\mathcal{S}^{m}_{i},\gamma_{i},\mathcal{S}_{-x}),$

a derivation of which can be found in Section B.2. To now predict $x$, one can
for example take the maximum a posteriori (MAP) estimate

$\hat{x}=\operatorname*{argmax}\limits_{x\in\mathcal{V}}\hat{p}(x|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n},\mathcal{S}_{-x}).$

In this simulation, we considered assessing the agreement of
$\hat{p}(x|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n},\mathcal{S}_{-x})$ with
$p(x|\mathcal{S}^{m},\mathcal{S}_{-x})$ as $n$ grows by comparing their
predictive accuracy on a common test sample. We adopted the same modelling
set-up as Section 5.1, jointly varying the dispersion and sample size, in this
case considering $\gamma_{\text{true}}=3.7,4.2,4.5,4.9$ and $n=25,50,75,100$.
However, in a slight deviation we here re-sampled the mode in each repetition
from a fixed Hollywood model.

Now, for a given pair $(\gamma_{\text{true}},n)$ and a pre-specified number of
test samples $n_{\text{test}}$, in a single repetition we (i) sampled mode
$\mathcal{S}_{\text{true}}\sim\text{Hollywood}(\alpha,-\alpha V,\nu)$, with
$\alpha=-0.35$ ($V$ and $\nu$ as in Sections 5.1 and 5.2) (ii) sampled
training and testing data $\\{\mathcal{S}^{(i)}\\}_{i=1}^{n+n_{\text{test}}}$
from an $\text{SIS}(\mathcal{S}_{\text{true}},\gamma_{\text{true}})$ model,
(iii) obtained a sample $\\{(\mathcal{S}^{m}_{i},\gamma_{i})\\}_{i=1}^{m}$
from the posterior
$p(\mathcal{S}^{m},\gamma|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})$, that is, using
the $n$ training samples, (iv) for each $i=n+1,\dots,n+n_{\text{test}}$ and
for each entry of $\mathcal{S}^{(i)}$ (that is, each entry of each
interaction) we assumed it to be missing and predicted its value via the MAP
estimate of both
$\hat{p}(x|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n},\mathcal{S}_{-x})$ and
$p(x|\mathcal{S}^{m},\mathcal{S}_{-x})$, and finally (v) returned the
proportion of times each prediction was correct. For (ii) and (iii) we used
the same MCMC schemes as previous simulations.

Figure 10 summarises the output of 100 repetitions for each pair
$(\gamma_{\text{true}},n)$, with $n_{\text{test}}=100$ in each repetition. For
each $\gamma_{\text{true}}$, we see the accuracy of the posterior predictive
is typically lower than the true predictive when the number of training
samples is small ($n=25$), with this difference in accuracy diminishing as $n$
grows, so that by $n=100$ training samples the accuracy of the posterior
predictive is almost indistinguishable from that of the true predictive. With
similarity in predictive accuracy serving as a proxy for the coherence of the
posterior predictive and the true data generating distribution, it thus
appears the posterior predictive is more closely resembling the true data
generating distribution as the number of training samples grows, as desired.
One can also observe influence from variability as controlled by
$\gamma_{\text{true}}$, where for larger $\gamma_{\text{true}}$ the predictive
accuracy is often higher (for both the posterior and the true predictive), as
one might expect since the test data in such cases will look more like the
true mode and thus more entries thereof will be easily predicted.

Figure 10: Summary of posterior predictive simulation (Section 5.3). In each
subplot, overlaid markers show the observed predictive accuracy being
summarised by boxplots, with the four leftmost boxes concerning the posterior
predictive at different sample sizes, whilst the rightmost box shows the
accuracy of the true predictive across all sample sizes (collated together).

##  6 Real Data Analysis

In this section, we illustrate the applicability of our method by analysing
the Foursquare data set of Yang et al. (2015). As mentioned in Section 1, an
alternative approach to ours is to first aggregate observations to form graphs
before applying a suitable graph-based method. As such, we compare our
inference with some graph-based estimates. Note that in aggregating
observations to form graphs one implicitly makes the assumption that the order
of interaction arrival is irrelevant. Hence, for fairness we opt to make this
comparison with our SIM model.

###  6.1 Data Background and Processing

For this analysis we looked at one month of check-in data (over the period
from 12 April to 12 May 2012, from the New York and Tokyo data
set222https://sites.google.com/site/yangdingqi/home/foursquare-dataset.),
focusing in particular on those in New York. Each check-in event consists of a
(i) user id, identifying which user enacted the check-in (anonymised) (ii)
venue id, unique to each venue, (iii) venue category, (iv) latitutde and
longitude, and (v) timestamp.

As discussed in Section 1, we view this as interaction network data by seeing
a day of check-ins for a single user as a path through the venue categories.
Note the venue category labels have a hierarchical structure, with those given
by Yang et al. (2015) being the low-level. For example, the category “Jazz
Club” is a subcategory of “Music Venue”, which is itself a subcategory of
“Arts & Entertainment”. In this analysis, we opted to use the highest-level
venue categories, so for example here we would use “Arts & Entertainment”.

Before proceeding with our analysis, we further filtered the data. Firstly, it
is clearly possible a user might only check-in to a single venue on a given
day. Since our analysis is based on interaction multisets and concerns the
movements of users between venue categories, such observations provide little
information. Furthermore, such observations will be disregarded when
aggregating to form graphs, and therefore would not feature in any of the
graph-based approaches with which we intend to compare. As such, we considered
only days where a user had checked-in to at least two venues. To further
ensure each observation contained enough information, we considered only users
with at least 10 observed days of check-ins. This left a total of 402
observations, from which we extracted a subset of 100 to analyse using a
criterion based upon the distance metric used in our model fit (further
details in Section S5.1).

###  6.2 SIM Model Fit

Following data processing, we were left with a sample of multisets
$\\{\mathcal{E}^{(i)}\\}_{i=1}^{n}$, where each
$\mathcal{E}^{(i)}=\left\\{\mathcal{I}^{(i)}_{1},\dots,\mathcal{I}^{(i)}_{N^{(i)}}\right\\}$
denotes the data of the $i$th user, with $\mathcal{I}^{(i)}_{j}$ denoting a
single day of their check-ins. Recalling the inferential questions of interest
outlined in Section 1, we now use our methodology to obtain (a) an average
multiset of paths, and (b) a measure of variability.

In particular, using the Bayesian inference approach outlined in Appendix S4,
we fit our SIM model to these data. We made use of the matching distance
$d_{\text{M}}$, with the LSP distance $d_{\text{LSP}}$ between paths. As
discussed in Section 3.5, a consequence of assuming this distance is that our
inferred mode will contain paths often appearing as subpaths in the observed
data. For our priors, we assumed
$\mathcal{E}^{m}\sim\text{SIM}(\hat{\mathcal{E}},3.0)$, with
$\hat{\mathcal{E}}$ denoting the sample Frechét mean of the observed data
$\\{\mathcal{E}^{(i)}\\}_{i=1}^{n}$, whilst we assumed
$\gamma\sim\text{Gamma}(5,1.67)$. Via our MCMC scheme, we then obtained a
sample $\\{(\mathcal{E}_{i}^{m},\gamma_{i})\\}_{i=1}^{M}$, from the posterior
$p(\mathcal{E}^{m},\gamma|\\{\mathcal{E}^{(i)}\\}_{i=1}^{n})$, obtaining a
total of $M=500$ samples with a burn-in period of 25,000 and taking a lag of
50 between samples.

Given the posterior sample $\\{(\mathcal{E}_{i}^{m},\gamma_{i})\\}_{i=1}^{M}$,
we subsequently obtained point estimates
$(\hat{\mathcal{E}}^{m},\hat{\gamma})$, with the mode estimate
$\hat{\mathcal{E}}^{m}$ functioning as our desired average, and $\hat{\gamma}$
a measure of data variability. In particular, we considered the following

$\displaystyle\hat{\mathcal{E}}^{m}=\operatorname*{argmin}_{\mathcal{E}\in\\{\mathcal{E}^{m}_{i}\\}_{i=1}^{M}}\sum_{i=1}^{M}d^{2}_{\text{M}}(\mathcal{E}^{m}_{i},\mathcal{E})$
$\displaystyle\hat{\gamma}=\frac{1}{M}\sum_{i=1}^{M}\gamma_{i}$

that is, the Frechét mean for the mode and the arithmetic mean for the
dispersion, both obtained from their respective posterior samples.

As mentioned, due to our choice of distance, the inferred mode
$\hat{\mathcal{E}}^{m}$ represents a collection of pathways frequently seen
together in the observed data. To visualise this, we consider plotting the
paths of $\hat{\mathcal{E}}^{m}$ alongside those of its two nearest
observations. Supposing the data points have been labelled such that
$\mathcal{E}^{(1)}$ and $\mathcal{E}^{(2)}$ denote the first and second
nearest neighbours of $\hat{\mathcal{E}}^{m}$ with respect to $d_{\text{M}}$,
writing these as follows

$\displaystyle\hat{\mathcal{E}}^{m}=\left\\{\hat{\mathcal{I}}^{m}_{1},\dots,\hat{\mathcal{I}}^{m}_{N^{m}}\right\\}$
$\displaystyle\mathcal{E}^{(1)}=\left\\{\mathcal{I}^{(1)}_{1},\dots,\mathcal{I}^{(1)}_{N_{1}}\right\\}$
$\displaystyle\mathcal{E}^{(2)}=\left\\{\mathcal{I}^{(2)}_{1},\dots,\mathcal{I}^{(2)}_{N_{2}}\right\\},$

Figures 11, 12 and 13 visualise the paths of $\hat{\mathcal{E}}^{m}$,
$\mathcal{E}^{(1)}$ and $\mathcal{E}^{(2)}$ alongside one another. In each,
the paths of $\mathcal{E}^{(1)}$ and $\mathcal{E}^{(2)}$ have been aligned in
accordance with the optimal matching found when evaluating their distance from
$\hat{\mathcal{E}}^{m}$ via $d_{\text{M}}$. In particular, in the $j$th row we
plot $\hat{\mathcal{I}}^{m}_{j}$ alongside $\mathcal{I}^{(1)}_{j}$ and
$\mathcal{I}^{(2)}_{j}$, denoting the paths matched to
$\hat{\mathcal{I}}^{m}_{j}$ when evaluating
$d_{\text{M}}(\hat{\mathcal{E}}^{m},\mathcal{E}^{(1)})$ and
$d_{\text{M}}(\hat{\mathcal{E}}^{m},\mathcal{E}^{(2)})$, respectively. The
paths of $\mathcal{E}^{(1)}$ and $\mathcal{E}^{(2)}$ not matched to any of
$\hat{\mathcal{E}}^{m}$ are then shown in the remaining rows.

Here one can observe paths of $\hat{\mathcal{E}}^{m}$ do indeed appear as
subpaths within those of $\mathcal{E}^{(1)}$ and $\mathcal{E}^{(2)}$. In fact,
in first three rows of Figure 11 they are equivalent, that is
$\hat{\mathcal{I}}^{m}_{j}=\mathcal{I}^{(1)}_{j}=\mathcal{I}^{(2)}_{j}$,
whilst for the remaining rows of Figure 11 and those of Figures 12 and 13 we
begin to see differences in the observed paths relative to those of the
estimated mode, however, in almost all cases, the paths in the mode
$\hat{\mathcal{I}}^{m}_{j}$ continue to feature as subpaths of both
$\mathcal{I}^{(1)}_{j}$ and $\mathcal{I}^{(2)}_{j}$. Note also that no paths
in $\hat{\mathcal{E}}^{m}$ are of length greater than two. At face value this
might seem to imply use of this method gains nothing over a graph-based
approach. However, as we illustrate in the next section, the subtle difference
is that our inference is unambiguous.

For the dispersion, we have $\hat{\gamma}\approx 3.02$, with a trace-plot of
the posterior samples $\\{\gamma_{i}\\}_{i=1}^{M}$ from which this estimate
was obtained shown in the left-hand plot of Figure 14. To aide interpretation
of $\hat{\gamma}$, the right hand plot of Figure 14 visualises the
distribution of $d_{\text{M}}(\mathcal{E},\hat{\mathcal{E}}^{m})$ where
$\mathcal{E}\sim\text{SIM}(\hat{\mathcal{E}}^{m},\gamma)$ for different values
of $\gamma$, each boxplot summarising 1,000 samples drawn from the respective
multiset model via our iMCMC algorithm (Section S4.5). A comparison with our
estimate $\hat{\gamma}$ shows that we expect the distance of samples to the
mode to be around 25 (from $\gamma=3.0$), which, since we used the matching
distance, can be seen as the average number of edit operations required to
transform the mode into an observation. Considering the mode has 18 entries in
total (9 paths of length two), this implies a reasonable amount of variability
in the observed data.

Figure 11: A subset of paths from our point estimate $\hat{\mathcal{E}}^{m}$
for the Foursquare data, alongside those of $\mathcal{E}^{(1)}$ and
$\mathcal{E}^{(2)}$, its two nearest neighbours. Paths are aligned according
to the optimal matching found when evaluating
$d_{\text{M}}(\hat{\mathcal{E}}^{m},\mathcal{E}^{(i)})$ for each neighbour
$\mathcal{E}^{(i)}$. For each observed path $\mathcal{I}^{(i)}_{j}$, dashed
pink edges and pink vertices indicate differences with
$\hat{\mathcal{I}}^{m}_{j}$, with edges labels indicating the order of vertex
visits. The remaining paths can be seen in Figures 12 and 13. Figure 12: Paths
of our point estimate $\hat{\mathcal{E}}^{m}$ for the Foursquare data,
alongside those of $\mathcal{E}^{(1)}$ and $\mathcal{E}^{(2)}$, its two
nearest neighbours. The remaining paths can be seen in Figures 11 and 13.
Figure 13: Paths of our point estimate $\hat{\mathcal{E}}^{m}$ for the
Foursquare data, alongside those of $\mathcal{E}^{(1)}$ and
$\mathcal{E}^{(2)}$, its two nearest neighbours. The remaining paths can be
seen in Figures 11 and 12. Figure 14: Summary of inference for the dispersion
for the Foursquare data. Left shows a trace-plot of the posterior samples
$\\{\gamma_{i}\\}_{i=1}^{m}$, whilst the right plot summarises the
distribution of distances to the inferred mode for different values of
$\gamma$, aiding interpretation of our estimate $\hat{\gamma}$.

###  6.3 Comparison with Graph-Based Inferences

An alternative to applying our methodology is to convert observations to
graphs in a pre-processing step, before applying a suitable graph-based
method. As such, we consider striking a comparison between this approach and
ours. The intention here is twofold. On one hand, to show the graph-based
inferences are not too dissimilar from that obtained via our approach. Whilst
on the other, that our approach goes beyond the graph-based methods, in so far
as producing an inference which is unambiguous regarding the presence of
higher-order information in the observed data.

Given the observed sample $\\{\mathcal{E}^{(i)}\\}_{i=1}^{n}$, one can obtain
a sample of graphs $\\{\mathcal{G}^{(i)}\\}_{i=1}^{n}$ via aggregation,
namely, by letting $\mathcal{G}^{(i)}=\mathcal{G}_{\mathcal{E}^{(i)}}$, the
graph obtained by aggregating the paths of $\mathcal{E}^{(i)}$, as outlined in
Section 2. In the same way that our estimate $\hat{\mathcal{E}}^{m}$
summarises the sample $\\{\mathcal{E}^{(i)}\\}_{i=1}^{n}$, one can now
consider obtaining a graph $\hat{\mathcal{G}}$ which summarises the sample
$\\{\mathcal{G}^{(i)}\\}_{i=1}^{n}$. This can be achieved through a variety of
different approaches, the choice of which will depend on whether the
$\mathcal{G}^{(i)}$ are graphs or multigraphs. We will consider both cases
here. In each instance, the ouput summary $\hat{\mathcal{G}}$ will be either a
graph or multigraph, which we then compare with
$\mathcal{G}_{\hat{\mathcal{E}}^{m}}$, the aggregation of our estimate
$\hat{\mathcal{E}}^{m}$.

To aide this exposition, we will make use of the graph adjacency matrix. For a
graph $\mathcal{G}=(\mathcal{E},\mathcal{V})$, where
$\mathcal{V}=\\{1,\dots,V\\}$, its adjacency matrix
$A^{\mathcal{G}}\in\\{0,1\\}^{V\times V}$ is the binary matrix with

$\displaystyle A^{\mathcal{G}}_{ij}=\begin{cases}1&\text{if
}(i,j)\in\mathcal{E}\\\ 0&\text{else},\end{cases}$

whilst if $\mathcal{G}=(\mathcal{E},\mathcal{V})$ is a multigraph its
adjacency matrix $A^{\mathcal{G}}\in\mathbb{Z}_{\geq 0}^{V\times V}$, is
defined by letting $A^{\mathcal{G}}_{ij}$ equal the number of times $(i,j)$
appears in $\mathcal{E}$. Note there is a one-to-one correspondence between
graphs and adjacency matrices, and as such they can be used interchangeably
for convenience.

In the case where each $\mathcal{G}^{(i)}$ is a graph, and thus each
$A^{\mathcal{G}^{(i)}}$ is a binary matrix, a simple model-free summary is the
majority vote, which we denote $\hat{\mathcal{G}}_{\text{MV}}$, where we
include an edge if it was observed in at least one half of the observations.
More formally, $\hat{\mathcal{G}}_{\text{MV}}$ can be defined in terms of its
adjacency matrix as follows

$A^{\hat{\mathcal{G}}_{\text{MV}}}_{ij}=\mathds{1}(\bar{A}_{ij}\geq 0.5),$

where $\bar{A}$ is the real-valued matrix with entries
$\bar{A}_{ij}=\frac{1}{n}\sum_{k=1}^{n}A_{ij}^{\mathcal{G}^{(k)}}$, that is,
the entry-wise average of the observed adjacency matrices. As a model-based
alternative, we turn to the centered Erdös-Rényi (CER) model of Lunagómez et
al. (2021). Using the notation
$\mathcal{G}\sim\text{CER}(\mathcal{G}^{m},\alpha)$ when a graph $\mathcal{G}$
was drawn from the CER model with mode $\mathcal{G}^{m}$ (a graph), and noise
parameter $0\leq\alpha\leq 0.5$, we assumed the following hierachical model

$\displaystyle\mathcal{G}^{(i)}\,|\,\mathcal{G}^{m},\alpha$
$\displaystyle\sim\text{CER}(\mathcal{G}^{m},\alpha)\quad\text{(for
$i=1,\dots,n$)}$ $\displaystyle\mathcal{G}^{m}$
$\displaystyle\sim\text{CER}(\mathcal{G}_{0},\alpha_{0})$
$\displaystyle\alpha$ $\displaystyle\sim
0.5\cdot\text{Beta}(\beta_{1},\beta_{2})$

where $\mathcal{G}_{0}$ (a graph), $0\leq\alpha_{0}\leq 0.5$, $\beta_{1}>0$
and $\beta_{2}>0$ denote hyperparameters. For this analysis, we assumed an
uniformative set-up, letting $\mathcal{G}_{0}=\hat{\mathcal{G}}_{\text{MV}}$
and $\alpha_{0}=0.5$, leading to a uniform distribution over the space of
graphs for the prior on $\mathcal{G}^{m}$, whilst we took
$\beta_{1}=\beta_{2}=1$, similarly leading to the uniform distribution over
the interval $(0,0.5)$ for the prior on $\alpha$. Following the scheme of
Lunagómez et al. (2021), we drew a sample
$\\{(\mathcal{G}^{m}_{i},\alpha_{i})\\}_{i=1}^{M}$ from the posterior
$p(\mathcal{G}^{m},\alpha|\\{\mathcal{G}^{(i)}\\}_{i=1}^{n})$ via MCMC,
obtaining the desired summary via the sample Frechét mean

$\hat{\mathcal{G}}_{\text{CER}}=\operatorname*{argmin}_{\mathcal{G}\in\\{\mathcal{G}^{m}_{i}\\}}\sum_{i=1}^{n}d^{2}_{\text{H}}(\mathcal{G},\mathcal{G}^{m}_{i})$

where $d_{\text{H}}$ denotes the Hamming distance between graphs (Lunagómez et
al., 2021; Donnat and Holmes, 2018). Figures 15(a) and 15(b) show these two
un-weighted summaries, $\hat{\mathcal{G}}_{\text{CER}}$ and
$\hat{\mathcal{G}}_{\text{MV}}$, respectively, for the Foursquare data, where
it transpires that
$\hat{\mathcal{G}}_{\text{CER}}=\hat{\mathcal{G}}_{\text{MV}}$.

In the case where each $\mathcal{G}^{(i)}$ is a multigraph, and thus each
$A^{\mathcal{G}^{(i)}}$ is a matrix of non-negative integers, an analogous
model-free summary can be obtained by rounding the entries of $\bar{A}$ to the
nearest integer. Referring to this as the rounded mean estimate and denoting
it $\hat{\mathcal{G}}_{\text{RM}}$, it can be defined formally via its
adjacency matrix as follows

$A_{ij}^{\hat{\mathcal{G}}_{\text{RM}}}=\lfloor\bar{A}_{ij}\rfloor+\mathds{1}(\bar{A}_{ij}-\lfloor\bar{A}_{ij}\rfloor\geq
0.5),$

where the notation $\lfloor x\rfloor$ for $x\in\mathbb{R}$ denotes the floor
function. As a model-based approach, we consider using the SNF models proposed
by Lunagómez et al. (2021). Though originally proposed to model graphs, they
can be readily extended to handle multigraphs (see Section S5.2). Use of the
SNF, like our models, requires specification of a distance metric between
graphs. We considered taking the absolute difference of edge multiplicities,
or alternatively, the adjacency matrix entries, that is

$d_{1}(\mathcal{G},\mathcal{G}^{\prime})=\sum_{i,j}\left|A^{\mathcal{G}}_{ij}-A^{\mathcal{G}^{\prime}}_{ij}\right|,$

which can be seen as the generalisation of the Hamming distance to
multigraphs. Adopting the notation
$\mathcal{G}\sim\text{SNF}(\mathcal{G}^{m},\gamma)$ when a graph $\mathcal{G}$
is drawn from the SNF model with mode $\mathcal{G}^{m}$ (a multigraph) and
dispersion $\gamma>0$, we assumed the following hierarchical model

$\displaystyle\mathcal{G}^{(i)}\,|\,\mathcal{G}^{m},\gamma$
$\displaystyle\sim\text{SNF}(\mathcal{G}^{m},\gamma)\quad\text{(for
$i=1,\dots,n$)}$ $\displaystyle\mathcal{G}^{m}$
$\displaystyle\sim\text{SNF}(\mathcal{G}_{0},\gamma_{0})$
$\displaystyle\gamma$ $\displaystyle\sim\text{Gamma}(\alpha,\beta)$

where $\mathcal{G}_{0}$ (a multigraph), $\gamma_{0}>0$, $\alpha>0$ and
$\beta>0$ are hyperparameters. For this analysis, we took $\mathcal{G}_{0}$ to
be the sample Frechét mean of the observed multigraphs
$\\{\mathcal{G}^{(i)}\\}_{i=1}^{n}$ with respect to the distance $d_{1}$,
whilst we let $\gamma_{0}=0.1$, $\alpha=3$ and $\beta=1$. Again, we obtained a
sample $\\{(\mathcal{G}^{m}_{i},\gamma_{i})\\}_{i=1}^{M}$ from the posterior
$p(\mathcal{G}^{m},\gamma|\\{\mathcal{G}^{(i)}\\}_{i=1}^{n})$ via MCMC, before
invoking the sample Frechét mean to obtain the desired summary

$\hat{\mathcal{G}}_{\text{SNF}}=\operatorname*{argmin}_{\mathcal{G}\in\\{\mathcal{G}^{(i)}\\}_{i=1}^{n}}\sum_{i=1}^{n}d_{1}^{2}(\mathcal{G},\mathcal{G}^{m}_{i}).$

Note that the posterior here will be doubly-intractable, necessitating use of
a specialised MCMC algorithm. Lunagómez et al. (2021) adopted the algorithm of
Møller et al. (2006), however, since here we consider multigraphs, we cannot
apply their scheme directly. Instead, we took an alternative approach via the
exchange algorithm (Murray et al., 2006), details of which can be found in
Section S5.2. Visualisations of these two multigraph summaries,
$\hat{\mathcal{G}}_{\text{SNF}}$ and $\hat{\mathcal{G}}_{\text{RM}}$, can be
seen in Figures 15(c) and 15(d), respectively.

Comparing the graph-based methods amongst themselves, we see a slight
variation in the signal they uncover. For example, in taking edge
multiplicities into account, the multigraph-based estimate
$\hat{\mathcal{G}}_{\text{RM}}$ introduces edges which did not appear in
either of the graphs $\hat{\mathcal{G}}_{\text{CER}}$ and
$\hat{\mathcal{G}}_{\text{CER}}$, generally involving the node corresponding
to food venues. Conversely, the SNF-based estimate
$\hat{\mathcal{G}}_{\text{SNF}}$ appears to instead disregard edges which
appear in $\hat{\mathcal{G}}_{\text{CER}}$ and
$\hat{\mathcal{G}}_{\text{CER}}$. Sitting somewhere in between is the
aggregate of our estimate $\mathcal{G}_{\hat{\mathcal{E}}^{m}}$, which appears
to have some degree of similarity with each of the graph-based summaries, as
we intended to confirm. Moreover, a common theme seems to appear (for all
summaries). Namely, that visits to food venues feature strongly, often
followed or preceded by a visit to another food venue or some other venue
category, with shopping venues being a prevalent choice.

Naturally, one is inclined to ask if the aggregation of our estimate
$\mathcal{G}_{\hat{\mathcal{E}}^{m}}$ is not too dissimilar to those obtained
via analysis of the data processed into graphs or multigraphs, then what does
one gain by taking our approach? Recalling that we estimated
$\hat{\mathcal{E}}^{m}$, a multiset of paths, we argue this contains more
information regarding the signals present in the observed data than any graph-
based summary, assuming the data were truly path-observed. This comes back to
a point made in Section 2, namely that when one aggregates paths a loss of
information is incurred, which will in turn limit the conclusions one can draw
concerning the original data. For example, consider the CER-based summary
$\hat{\mathcal{G}}_{\text{CER}}$ of Figure 15(a), where we see the following
two edges

$\displaystyle e_{1}=(\text{recreation},\text{food})$ $\displaystyle
e_{2}=(\text{food},\text{shops}).$

This could imply at least two things. Perhaps many users went from a
recreation venue to a food venue, and separately, that is, on a different day,
from a food venue to a shopping venue. Alternatively, maybe many users traced
the path recreation $\to$ food $\to$ shops in a single day. Both are
possibilities, and from a graph-based summary there is no way of knowing which
is the case. However, in directly estimating a collection of paths, we can
make such distinctions. For example, considering our estimate
$\hat{\mathcal{E}}^{m}$ for the Foursquare data (Figures 11, 12 and 13), it
appears we are in the former case, since no paths therein are of length
greater than two.

(a) $\hat{\mathcal{G}}_{\text{CER}}$

(b) $\hat{\mathcal{G}}_{\text{MV}}$

(c) $\hat{\mathcal{G}}_{\text{SNF}}$

(d) $\hat{\mathcal{G}}_{\text{RM}}$

(e) $\mathcal{G}_{\hat{\mathcal{E}}^{m}}$

Figure 15: A comparison with graph-based inferences. Here (e) shows
$\mathcal{G}_{\hat{\mathcal{E}}^{m}}$, the aggregate multigraph of our point
estimate $\hat{\mathcal{E}}^{m}$ of Section 6.2, whilst (a)-(d) show
alternative inferences obtained via graph-based approaches outlined in Section
6.3. Note that (a) and (b) are graphs, whilst (c)-(e) are multigraphs, with
edge thickness proportional to their weight.

##  7 Discussion

In this paper, we have motivated and instantiated the study of multiple
interaction-network data. We have proposed a flexible Bayesian modelling
framework capable of analysing such data without the need to perform any
aggregation of observations. This has been supplemented with specialised MCMC
schemes, facilitating inference for our proposed models. Through simulation
studies, we have confirmed the efficacy of our approach and inference scheme,
whilst the applicability of our methodology has been illustrated by an
analysis of Foursquare check-in data, where we illustrated how our methodology
can be used to answer inferential questions (a) and (b) posed in Section 1.
Moreover, in comparing with graph-based methods we highlighted the extra
information one subtly gains by taking our approach.

Regarding future work, there are a few ways one might consider building upon
what has been proposed here. Firstly, a natural extension of our models is to
consider a mixture model, with our SIS or SIM models functioning as mixture
components, which would allow one to capture heterogeneity in the
observations, opening the door to answering question (c) of Section 1.
Secondly, on a more pragmatic note, one could also take steps to scale-up our
approach computationally. For example, one might be able to circumvent the
need to use the exchange algorithm if the normalising constant for a
particular distance metric was derived, as was the case for the CER model in
Lunagómez et al. (2021). Finally, if one is able to make an exchangeability
assumption for each observation, that is, the order in which paths arrive is
not of interest, then a slightly modified model structure could be considered,
reminiscent of the latent Dirichlet allocation (LDA) model (Blei et al.,
2003). Namely, one could assume each observation was drawn from some mixture
distribution over paths, with mixture components being shared between
observations but mixture proportions differing. This would also have a natural
non-parametric extension via the hierarchical Dirichlet process (HDP) (Teh et
al., 2006). It would be interesting to see how the inferences from such an
approach compare with ours, at least qualitatively, and whether any
computational benefit would be achieved.

More tangentially, one could also follow the path laid in the wider literature
on multiple networks and consider extending models designed to analyse a
single interaction network, for example, the models of Crane and Dempsey
(2018) or Williamson (2016).

Finally, it could be interesting to consider the situation where one has
access to covariate information at the level of observations. For example,
considering the Foursquare data, one might have additional information for
each user, such as their line of work or country of residence. Interest might
then be in defining a modelling framework which could be invoked to examine
for a relationship between covariates and observed data. Such developments
would mirror those in the wider literature on multiple-network data, such as
work on hypothesis testing (Ginestet et al., 2017; Durante and Dunson, 2018;
Ghoshdastidar et al., 2020; Chen et al., 2021).

##  Appendix A Sample Spaces

Here we provide extra details regarding the sample spaces for our SIS and SIM
models of Definition 1 and 2. Recall that distributions within the defined
families constitute distributions over sequences and multisets (of paths)
respectively. Thus, their respective sample spaces consist of all such
objects.

Assuming the vertex set $\mathcal{V}$ is fixed, we first define the space of
all paths over $\mathcal{V}$ as follows

$\displaystyle\bm{\mathcal{I}}^{*}:=\\{(x_{1},\dots,x_{n})\,:\,x_{i}\in\mathcal{V},n\geq
1\\},$

subsequently defining the space of interaction sequences
$\bm{\mathcal{S}}^{*}$ via

$\displaystyle\bm{\mathcal{S}}^{*}:=\\{(\mathcal{I}_{1},\dots,\mathcal{I}_{N})\,:\,\mathcal{I}_{i}\in\bm{\mathcal{I}}^{*},N\geq
1\\},$

whilst, with $\mathcal{E}_{\mathcal{S}}$ denoting the multiset obtained from
the sequence $\mathcal{S}$ by disregarding the order of paths therein, the
space of interaction multisets $\bm{\mathcal{E}^{*}}$ can be defined via

$\displaystyle\bm{\mathcal{E}^{*}}:=\\{\mathcal{E}_{\mathcal{S}}\,:\,\mathcal{S}\in\bm{\mathcal{S}}^{*}\\}$

where here we are abusing notation slightly, since we can have
$\mathcal{E}_{\mathcal{S}}=\mathcal{E}_{\mathcal{S}^{\prime}}$ for
$\mathcal{S}\not=\mathcal{S}^{\prime}$ (when equal up to a permutation of
interactions), but we just assume such values have been included once and so
$\bm{\mathcal{E}^{*}}$ is indeed a set.

We note that $\bm{\mathcal{E}^{*}}$ also admits another interpretation as a
partitioning of $\bm{\mathcal{S}}^{*}$ into equivalence classes. To see this,
first define an equivalence relation on $\bm{\mathcal{S}}^{*}$ via
permutations, in particular we write
$\mathcal{S}\overset{p}{\sim}\mathcal{S}^{\prime}$ if there is some
permutation $\sigma$ such that $\mathcal{S}^{\prime}=\mathcal{S}^{\sigma}$,
where
$\mathcal{S}^{\sigma}=(\mathcal{I}_{\sigma(1)},\dots,\mathcal{I}_{\sigma(N)})$
is the interaction sequence obtained by permuting the interactions of
$\mathcal{S}$ via $\sigma$. Now, observe that each
$\mathcal{E}\in\bm{\mathcal{E}^{*}}$ can be seen as an equivalence class of
interaction sequences obtained via $\overset{p}{\sim}$, that is

$\mathcal{E}=\\{\mathcal{S}\in\bm{\mathcal{S}}^{*}\,:\,\mathcal{S}\overset{p}{\sim}\tilde{\mathcal{S}}\\}$

where $\tilde{\mathcal{S}}$ denotes some arbitrary ordering of the
interactions of $\mathcal{E}$. Thus, $\bm{\mathcal{E}^{*}}$ is in a sense the
union of such sets and partitions $\bm{\mathcal{S}}^{*}$.

###  A.1 Bounding Dimensions

As highlighted Section 3, it can in practice be somewhat easier to work with
bounded sample spaces, since in the unbounded case our models are not
guaranteed to be proper, which can lead to a divergence of dimension when
doing MCMC sampling for some parameterisations (typically when $\gamma$ is
low). This can further cause computational instabilities if one uses such MCMC
sampling within other algorithms, as we do within our scheme to sample from
the posterior. In this section, we state our notation for the bounded sample
spaces.

With regards to the objects we consider, there are two things we can bound:
(i) the size of paths and (ii) the number of paths. Referring to these as the
inner and outer dimensions respectively, we specify two integers $K$ and $L$
bounding their values and define our sample spaces accordingly. Assuming that
the vertex set $\mathcal{V}$ is fixed, and $K\in\mathbb{Z}_{\geq 1}$ we let

$\displaystyle\bm{\mathcal{I}}^{*}_{K}:=\\{(x_{1},\dots,x_{n})\,:\,x_{i}\in\mathcal{V},1\leq
n\leq K\\}$

denote the space of paths up to length $K$, and then with
$L\in\mathbb{Z}_{\geq 1}$ we let

$\displaystyle\bm{\mathcal{S}}^{*}_{K,L}:=\\{(\mathcal{I}_{1},\dots,\mathcal{I}_{N})\,:\,\mathcal{I}_{i}\in\bm{\mathcal{I}}^{*}_{K},1\leq
N\leq L\\},$

denote the space of interaction sequences with at most $L$ paths of length at
most $K$. The analogous bounded space of interaction multisets in then given
by

$\displaystyle\bm{\mathcal{E}^{*}}_{K,L}:=\\{\mathcal{E}_{\mathcal{S}}\,:\,\mathcal{S}\in\bm{\mathcal{S}}^{*}_{K,L}\\}.$

##  Appendix B Simulation Studies

This section contains supporting details for the simulation studies of Section
5. In particular, we discuss how parameters were chosen for the simulation of
Section 5.2, and provide a derivation for the posterior predictive
approximation used in Section 5.3.

###  B.1 Posterior Concentration Parameter Choices

Recall that in the simulation of Section 5.2 we re-sampled the true mode via

$\mathcal{S}_{\text{true}}\sim\text{Hollywood}(\alpha,-\alpha V,\nu)$

where $V=20$ and $\nu=\text{TrPoisson}(3,1,10)$, whilst $\alpha<0$ we varied.
As mentioned in Section 5.2, the parameter $\alpha$ can be seen to control the
tail of the vertex count distribution. As such, rather than choosing $\alpha$
on an even grid we instead consult a summary measure quantifying the ‘heavy-
tailedness’ of the degree distribution, before choosing values so as to evenly
represent different structures for $\mathcal{S}_{\text{true}}$ (as quantified
by this degree distribution).

For a given observation $\mathcal{S}$, recall the following definition

$k_{\mathcal{S}}(v):=\text{\\# times $v$ appears in $\mathcal{S}$},$

which for each $\mathcal{S}$ implies a sample
$\\{k_{\mathcal{S}}(v)\,:v\in\mathcal{V},\,k_{\mathcal{S}}(v)>0\\}$, similar
to the degree distribution. Now, the summary measure we considered was the 95%
quantile of this sample.

Through simulation, we examined how $\alpha$ controls the expected value of
this 95% quantile (expected since $\mathcal{S}$ is sample randomly from a
Hollywood model). In particular, for a range of $\alpha$ values, we (i) drew a
sample $\\{\mathcal{S}^{(i)}\\}_{i=1}^{n}$, from a
$\text{Hollywood}(\alpha,-\alpha V,\nu)$ model, taking $\nu$ and $V$ as above,
drawing a total of $N=10$ paths in each case, then (ii) for $i=1,\dots,n$ we
evaluated the 95% quantile of the sample
$\\{k_{\mathcal{S}^{(i)}}(v)\,:v\in\mathcal{V},\,k_{\mathcal{S}^{(i)}}(v)>0\\}$,
before returning the mean value of these quantiles.

Figure 16 summarises the output with $n=1000$ samples, where circular markers
show the mean quantiles. Towards choosing simulation parameters, we next
constructed a function mapping all $\alpha<0$ to an expected quantile via a
linear interpolation, as shown in Figure 16 by the dashed line, which allowed
us to select $\alpha$ (red crosses) providing an even spread of expected
degree-distribution 95% qauntiles.

Figure 16: Summary of Hollywood model simulation used to select parameters for
simulation of Section 5.2. Plot shows simulated mean degree distribution 95%
quantiles for $\text{Hollywood}(\alpha,-\alpha V,\nu)$ model, where $V=20$,
$\nu=\text{TrPoisson}(3,1,10)$ and $\alpha$ varies. Via linear interpolation
(dashed line), we choose $\alpha$ values (crosses) to get an even spread over
the expected degree distribution quantiles.

###  B.2 Posterior Predictive for Missing Entries

Here we show how one can obtain an approximation for the missing-entry
posterior predictive using a sample from the posterior, as used in Section
5.3. First, observe that any sample
$\\{(\mathcal{S}^{m}_{i},\gamma_{i})\\}_{i=1}^{m}$ from the posterior implies
the following atomic approximation thereof

$\displaystyle\hat{p}(\mathcal{S}^{m},\gamma|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})$
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}\mathds{1}(\mathcal{S}^{m}=\mathcal{S}^{m}_{i})\cdot\delta(\gamma-\gamma_{i})$
(14)

where $\delta(\cdot)$ is the Dirac delta function.

As in Section 5.3, with $\mathcal{S}_{x}$ denoting the observation with
missing entry filled in to be $x$, then given some parameters
$(\mathcal{S}^{m},\gamma)$ we have the true predictive for $x$ given by

$\displaystyle
p(x|\mathcal{S}^{m},\gamma,\mathcal{S}_{-x}):=\frac{1}{Z(\mathcal{S}^{m},\gamma,\mathcal{S}_{-x})}\exp\\{-\gamma
d_{S}(\mathcal{S}_{x},\mathcal{S}^{m})\\}$

with

$Z(\mathcal{S}^{m},\gamma,\mathcal{S}_{-x})=\sum_{x\in\mathcal{V}}\exp\\{-\gamma
d_{S}(\mathcal{S}_{x},\mathcal{S}^{m})\\}.$

The posterior predictive is now obtained by averaging with respect to the
posterior

$\displaystyle p(x|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n},\mathcal{S}_{-x})$
$\displaystyle=\sum_{\mathcal{S}^{m}\in\bm{\mathcal{S}}^{*}}\int_{\mathbb{R}_{+}}p(x|\mathcal{S}^{m},\gamma,\mathcal{S}_{-x})p(\mathcal{S}^{m},\gamma|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})d\gamma,$

which we can now approximate by substituting in (14) as follows

$\displaystyle\hat{p}(x|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n},\mathcal{S}_{-x})$
$\displaystyle:=\sum_{\mathcal{S}^{m}\in\bm{\mathcal{S}}^{*}}\int_{\mathbb{R}_{+}}p(x|\mathcal{S}^{m},\gamma,\mathcal{S}_{-x})\hat{p}(\mathcal{S}^{m},\gamma|\\{\mathcal{S}^{(i)}\\}_{i=1}^{n})d\gamma$
$\displaystyle=\sum_{\mathcal{S}^{m}\in\bm{\mathcal{S}}^{*}}\int_{\mathbb{R}_{+}}p(x|\mathcal{S}^{m},\gamma)\left(\frac{1}{m}\sum_{i=1}^{m}\mathds{1}\left(\mathcal{S}^{m}=\mathcal{S}^{m}_{i}\right)\delta(\gamma-\gamma_{i})\right)d\gamma$
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}p(x|\mathcal{S}^{m}_{i},\gamma_{i}),$

which is exactly as stated in Section 5.3.

A pragmatic note here is that as the posterior concentrates the number of
unique values in the sample $\\{\mathcal{S}^{m}_{i}\\}_{i=1}^{m}$ will
typically not be too large. Since we need only evaluate the distance metric
(which is typically quite costly) at these values, this predictive is feasible
to evaluate.

## References

  * Arroyo et al. (2021) Jesús Arroyo, Avanti Athreya, Joshua Cape, Guodong Chen, Carey E Priebe, and Joshua T Vogelstein.  Inference for multiple heterogeneous networks with a common invariant subspace.  _Journal of machine learning research_ , 22(142), 2021.
  * Behrens and Sporns (2012) Timothy E.J. Behrens and Olaf Sporns.  Human connectomics.  _Current Opinion in Neurobiology_ , 22:144–153, 2012.
  * Blei et al. (2003) David M Blei, Andrew Y Ng, and Michael I. Jordan.  Latent dirichlet allocation.  _Journal of Machine Learning Research_ , 3:993–1022, 2003\.
  * Bolt et al. (2022) George Bolt, Simón Lunagómez, and Christopher Nemeth.  Distances for comparing multisets and sequences.  _arXiv preprint arXiv:2206.08858_ , 2022.
  * Cai et al. (2016) Diana Cai, Trevor Campbell, and Tamara Broderick.  Edge-exchangeable graphs and sparsity.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Caron and Fox (2017) François Caron and Emily B Fox.  Sparse graphs using exchangeable random measures.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 79(5):1295–1366, 2017.
  * Chen et al. (2021) Li Chen, Jie Zhou, and Lizhen Lin.  Hypothesis testing for populations of networks.  _Communications in Statistics-Theory and Methods_ , pages 1–24, 2021\.
  * Chung et al. (2021) Jaewon Chung, Eric Bridgeford, Jesús Arroyo, Benjamin D Pedigo, Ali Saad-Eldin, Vivek Gopalakrishnan, Liang Xiang, Carey E Priebe, and Joshua T Vogelstein.  Statistical connectomics.  _Annual Review of Statistics and Its Application_ , 8:463–492, 2021.
  * Crane and Dempsey (2018) Harry Crane and Walter Dempsey.  Edge exchangeable models for interaction networks.  _Journal of the American Statistical Association_ , 113:1311–1326, 2018.
  * Donnat and Holmes (2018) Claire Donnat and Susan Holmes.  Tracking network dynamics: A survey using graph distances.  _Annals of Applied Statistics_ , 12:971–1012, 2018.
  * Durante and Dunson (2018) Daniele Durante and David B Dunson.  Bayesian inference and testing of group differences in brain networks.  _Bayesian Analysis_ , 13(1):29–58, 2018.
  * Durante et al. (2017) Daniele Durante, David B Dunson, and Joshua T Vogelstein.  Nonparametric Bayes modeling of populations of networks.  _Journal of the American Statistical Association_ , 112:1516–1530, 2017.
  * Frank and Strauss (1986) Ove Frank and David Strauss.  Markov graphs.  _Journal of the American Statistical Association_ , 81(395):832–842, 1986.
  * Ghalebi et al. (2019a) Elahe Ghalebi, Hamidreza Mahyar, Radu Grosu, Graham W Taylor, and Sinead A Williamson.  A nonparametric Bayesian model for sparse dynamic multigraphs.  _arXiv preprint arXiv:1910.05098_ , 2019a.
  * Ghalebi et al. (2019b) Elahe Ghalebi, Hamidreza Mahyar, Radu Grosu, and Sinead Williamson.  Dynamic nonparametric edge-clustering model for time-evolving sparse networks.  _arXiv preprint arXiv:1905.11724_ , 2019b.
  * Ghoshdastidar et al. (2020) Debarghya Ghoshdastidar, Maurilio Gutzeit, Alexandra Carpentier, and Ulrike Von Luxburg.  Two-sample hypothesis testing for inhomogeneous random graphs.  _The Annals of Statistics_ , 48(4):2208–2229, 2020.
  * Ginestet et al. (2017) Cedric E. Ginestet, Jun Li, Prakash Balachandran, Steven Rosenberg, and Eric D. Kolaczyk.  Hypothesis testing for network data in functional neuroimaging.  _Annals of Applied Statistics_ , 11:725–750, 7 2017.
  * Gollini and Murphy (2016) Isabella Gollini and Thomas Brendan Murphy.  Joint modeling of multiple network views.  _Journal of Computational and Graphical Statistics_ , 25:246–265, 2016.
  * Hoff et al. (2002) Peter D Hoff, Adrian E Raftery, and Mark S Handcock.  Latent space approaches to social network analysis.  _Journal of the American Statistical Association_ , 97:1090–1098, 2002.
  * Holland and Leinhardt (1981) Paul W Holland and Samuel Leinhardt.  An exponential family of probability distributions for directed graphs.  _Journal of the American Statistical Association_ , 76(373):33–50, 1981.
  * Kuhn (1955) H. W. Kuhn.  The Hungarian method for the assignment problem.  _Naval Research Logistics Quarterly_ , 2:83–97, 1955.
  * Le et al. (2018) Can M Le, Keith Levin, and Elizaveta Levina.  Estimating a network from multiple noisy realizations.  _Electronic Journal of Statistics_ , 12:4697–4740, 2018\.
  * Lehmann and White (2021) Brieuc Lehmann and Simon White.  Bayesian exponential random graph models for populations of networks.  _arXiv preprint arXiv:2104.05110_ , 2021.
  * Levin et al. (2017) Keith Levin, Avanti Athreya, Minh Tang, Vince Lyzinski, and Carey E Priebe.  A central limit theorem for an omnibus embedding of multiple random dot product graphs.  In _2017 IEEE International Conference on Data Mining Workshops (ICDMW)_ , pages 964–967. IEEE, 2017.
  * Liang (2010) Faming Liang.  A double Metropolis-Hastings sampler for spatial models with intractable normalizing constants.  _Journal of Statistical Computation and Simulation_ , 80:1007–1022, 2010.
  * Lunagómez et al. (2021) Simón Lunagómez, Sofia C. Olhede, and Patrick J. Wolfe.  Modeling network populations via graph distances.  _Journal of the American Statistical Association_ , 116(536):2023–2040, 2021.
  * Mantziou et al. (2021) Anastasia Mantziou, Simon Lunagomez, and Robin Mitra.  Bayesian model-based clustering for multiple network data.  _arXiv preprint arXiv:2107.03431_ , 2021.
  * Mardia and Dryden (1999) K. V. Mardia and I. L. Dryden.  The complex Watson distribution and shape analysis.  _Journal of the Royal Statistical Society. Series B: Statistical Methodology_ , 61:913–926, 1999.
  * Marron and Dryden (2021) James Stephen Marron and Ian L. Dryden.  _Object Oriented Data Analysis (1st ed.)_.  Chapman and Hall/CRC, 2021.
  * Murray et al. (2006) Iain Murray, Zoubin Ghahramani, and David J.C. MacKay.  Mcmc for doubly-intractable distributions.  _Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence, UAI 2006_ , pages 359–366, 2006.
  * Møller et al. (2006) J. Møller, A N Pettitt, R Reeves, and K. K. Berthelsen.  An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants.  _Biometrika_ , 93:451–458, 2006.
  * Neklyudov et al. (2020) Kirill Neklyudov, Max Welling, Evgenii Egorov, and Dmitry Vetrov.  Involutive mcmc: A unifying framework.  In _International Conference on Machine Learning_ , pages 7273–7282. PMLR, 2020.
  * Newman (2018) Mark EJ Newman.  Estimating network structure from unreliable measurements.  _Physical Review E_ , 98(6):062321, 2018.
  * Nielsen and Witten (2018) Agnes Martine Nielsen and Daniela Witten.  The multiple random dot product graph model.  _arXiv preprint arXiv:1811.12172_ , 11 2018.
  * Nowicki and Snijders (2001) Krzysztof Nowicki and Tom A B Snijders.  Estimation and prediction for stochastic blockstructures.  _Journal of the American Statistical Association_ , 96(455):1077–1087, 2001.
  * Peixoto (2018) Tiago P Peixoto.  Reconstructing networks with unknown and heterogeneous errors.  _Physical Review X_ , 8(4):041011, 2018.
  * Peixoto and Rosvall (2017) Tiago P Peixoto and Martin Rosvall.  Modelling sequences and temporal networks with dynamic community structures.  _Nature Communications_ , 8(1):1–12, 2017.
  * Reyes and Rodriguez (2016) Perla Reyes and Abel Rodriguez.  Stochastic blockmodels for exchangeable collections of networks.  _arXiv preprint arXiv:1606.05277_ , 2016.
  * Scholtes (2017) Ingo Scholtes.  When is a network a network? multi-order graphical model selection in pathways and temporal networks.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1037–1046, 2017.
  * Stanley et al. (2016) Natalie Stanley, Saray Shai, Dane Taylor, and Peter J Mucha.  Clustering network layers with the strata multilayer stochastic block model.  _IEEE Transactions on Network Science and Engineering_ , 3(2):95–105, 2016.
  * Sweet et al. (2013) Tracy M Sweet, Andrew C Thomas, and Brian W Junker.  Hierarchical network models for education research: Hierarchical latent space models.  _Journal of Educational and Behavioral Statistics_ , 38(3):295–318, 2013.
  * Sweet et al. (2014) Tracy M Sweet, Andrew C Thomas, and Brian W Junker.  Hierarchical mixed membership stochastic blockmodels for multiple networks and experimental interventions.  _Handbook on Mixed Membership Models and their Applications_ , pages 463–488, 2014.
  * Teh et al. (2006) Yee Whye Teh, Michael I Jordan, Matthew J. Beal, and David M Blei.  Hierarchical dirichlet processes.  _Journal of the American Statistical Association_ , 101:1566–1581, 2006.
  * Veitch and Roy (2015) Victor Veitch and Daniel M Roy.  The class of random graphs arising from exchangeable random measures.  _arXiv preprint arXiv:1512.03099_ , 2015.
  * Vitelli et al. (2018) Valeria Vitelli, Øystein Sørensen, Marta Crispino, Arnoldo Frigessi, and Elja Arjas.  Probabilistic preference learning with the Mallows rank model.  _Journal of Machine Learning Research_ , 18:1–49, 2018\.
  * Wang et al. (2019) Shangsi Wang, Jesús Arroyo, Joshua T Vogelstein, and Carey E Priebe.  Joint embedding of graphs.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2019.
  * Williamson (2016) Sinead A Williamson.  Nonparametric network models for link prediction.  _The Journal of Machine Learning Research_ , 17(1):7102–7121, 2016.
  * Yang et al. (2015) Dingqi Yang, Daqing Zhang, Vincent W. Zheng, and Zhiyong Yu.  Modeling user activity preference by leveraging user spatial temporal characteristics in lbsns.  _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , 45:129–142, 2015.
  * Yin et al. (2022) Fan Yin, Weining Shen, and Carter T Butts.  Finite mixtures of ergms for modeling ensembles of networks.  _Bayesian Analysis_ , 1(1):1–39, 2022.
  * Young et al. (2022) Jean-Gabriel Young, Alec Kirkley, and MEJ Newman.  Clustering of heterogeneous populations of networks.  _Physical Review E_ , 105(1):014312, 2022.
  * Young and Scheinerman (2007) Stephen J Young and Edward R Scheinerman.  Random dot product graph models for social networks.  In _International Workshop on Algorithms and Models for the Web-Graph_ , pages 138–149. Springer, 2007.

## Supplementary Material

##  Appendix S1 Montonicity of the Entropy

Here we examine the entropy for the SIS and SIM model families (Definition 1
and 2 respectively), in particular, we confirm it is monotonic with respect to
the dispersion. A similar result was shown by Lunagómez et al. (2021)
(Proposition 3.3, proved pages 41-43 of Supplementary Material), and as such
we follow the reasoning used therein.

For the SIS model, recall the entropy is given by

$\displaystyle H(\mathcal{S}^{m},\gamma)$ $\displaystyle=-\mathbb{E}[\log
p(\mathcal{S}|\mathcal{S}^{m},\gamma)]$
$\displaystyle=-\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}\log\left(\frac{\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}}{Z(\mathcal{S}^{m},\gamma)}\right)\frac{\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}}{Z(\mathcal{S}^{m},\gamma)}$
$\displaystyle=-\bigg{(}\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\frac{\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}}{Z(\mathcal{S}^{m},\gamma)}$
$\displaystyle\qquad\qquad-\log
Z(\mathcal{S}^{m},\gamma)\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}\frac{\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}}{Z(\mathcal{S}^{m},\gamma)}\bigg{)}$
$\displaystyle=\gamma\left(\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}d_{S}(\mathcal{S},\mathcal{S}^{m})\frac{\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}}{Z(\mathcal{S}^{m},\gamma)}\right)+\log
Z(\mathcal{S}^{m},\gamma)$
$\displaystyle=\gamma\times\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})]+\log
Z(\mathcal{S}^{m},\gamma).$

Unfortunately, as was the case for the normalising constant
$Z(\mathcal{S}^{m},\gamma)$ (Section 3.3), since $\bm{\mathcal{S}}^{*}$ is
infinite we have no guarantee that $H(\mathcal{S}^{m},\gamma)$ will exist.
However, what we can say is that, when $H(\mathcal{S}^{m},\gamma)$ exists, it
is monotonic in $\gamma$. To show this, we first differentiate
$H(\mathcal{S}^{m},\gamma)$ with respect to $\gamma$

$\displaystyle\frac{\partial}{\partial\gamma}H(\mathcal{S}^{m},\gamma)=\frac{\partial}{\partial\gamma}\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})]+\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})]+\frac{\partial}{\partial\gamma}\log
Z(\mathcal{S}^{m},\gamma)$

where one has

$\displaystyle\frac{\partial}{\partial\gamma}\log Z(\mathcal{S}^{m},\gamma)$
$\displaystyle=\frac{\frac{\partial}{\partial\gamma}Z(\mathcal{S}^{m},\gamma)}{Z(\mathcal{S}^{m},\gamma)}$
(15)
$\displaystyle=\frac{1}{Z(\mathcal{S}^{m},\gamma)}\frac{\partial}{\partial\gamma}\left(\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}\exp\left\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\right\\}\right)$
$\displaystyle=\frac{1}{Z(\mathcal{S}^{m},\gamma)}\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}\frac{\partial}{\partial\gamma}\exp\left\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\right\\}$
$\displaystyle=\frac{1}{Z(\mathcal{S}^{m},\gamma)}\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}(-d_{S}(\mathcal{S},\mathcal{S}^{m}))\exp\left\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\right\\}$
$\displaystyle=-\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}d_{S}(\mathcal{S},\mathcal{S}^{m})\frac{1}{Z(\mathcal{S}^{m},\gamma)}\exp\left\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\right\\}$
$\displaystyle=-\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})],$

thus implying

$\displaystyle\frac{\partial}{\partial\gamma}H(\mathcal{S}^{m},\gamma)=\frac{\partial}{\partial\gamma}\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})].$

Now, we have

$\displaystyle\frac{\partial}{\partial\gamma}\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})]$
$\displaystyle=\frac{\partial}{\partial\gamma}\left(\frac{1}{Z(\mathcal{S}^{m},\gamma)}\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}d_{S}(\mathcal{S},\mathcal{S}^{m})\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}\right)$
$\displaystyle=-\frac{\frac{\partial}{\partial\gamma}Z(\mathcal{S}^{m},\gamma)}{Z(\mathcal{S}^{m},\gamma)^{2}}\left(\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}d_{S}(\mathcal{S},\mathcal{S}^{m})\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}\right)$
$\displaystyle\qquad\qquad-\frac{1}{Z(\mathcal{S}^{m},\gamma)}\left(\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}d_{S}(\mathcal{S},\mathcal{S}^{m})^{2}\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}\right)$
$\displaystyle=-\frac{\frac{\partial}{\partial\gamma}Z(\mathcal{S}^{m},\gamma)}{Z(\mathcal{S}^{m},\gamma)^{2}}\left(\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}d_{S}(\mathcal{S},\mathcal{S}^{m})\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}\right)$
$\displaystyle\qquad\qquad-\frac{1}{Z(\mathcal{S}^{m},\gamma)}\left(\sum_{\mathcal{S}\in\bm{\mathcal{S}}^{*}}d_{S}(\mathcal{S},\mathcal{S}^{m})^{2}\exp\\{-\gamma
d_{S}(\mathcal{S},\mathcal{S}^{m})\\}\right)$
$\displaystyle=-\frac{\frac{\partial}{\partial\gamma}Z(\mathcal{S}^{m},\gamma)}{Z(\mathcal{S}^{m},\gamma)}\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})]-\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})^{2}]$
(16)
$\displaystyle=\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})]^{2}-\mathbb{E}[d_{S}(\mathcal{S},\mathcal{S}^{m})^{2}]$
(17) $\displaystyle=\mathbb{V}\text{ar}[d_{S}(\mathcal{S},\mathcal{S}^{m})],$

where (17) follows from (16) by applying (15). Now, observe that if
$\frac{\partial}{\partial\gamma}H(\mathcal{S}^{m},\gamma)>0$ this implies
$H(\mathcal{S}^{m},\gamma)$ is monotonic in $\gamma$, as desired. By the
derivations above, this is equivalent to saying we have montonicity provided
$\mathbb{V}\text{ar}[d_{S}(\mathcal{S},\mathcal{S}^{m})]>0$. This result,
identical to that of Lunagómez et al. (2021), essentially says we have
monotonicity of the entropy with respect to $\gamma$ provided our distribution
is not a point mass.

Similar derivations can be obtained for the multiset models (Definition 2) by
a simple change of notation. For brevity, we do not repeat this here.

##  Appendix S2 The iExchange Algorithm

In this section, we outline the iExchange algorithm (Algorithm 1), a
generalisation of exchange algorithm (Murray et al., 2006) obtained by
incorporating the proposal generating mechanism of the iMCMC algorithm
(Neklyudov et al., 2020). As we show, the iExchange algorithm is itself an
iMCMC algorithm (with a particular form of involution), providing the
necessary theoretical justification. For completeness, we give background
details regarding both the exchange and iMCMC algorithms, before showing how
they can be combined.

Input: target density $p(\theta|\mathbf{x})\propto
p(\theta)\gamma(\mathbf{x}|\theta)/Z(\theta)$

Input: auxiliary density $q(u|\theta)$

Input: involution $f(\theta,u)$, i.e. $f^{-1}(\theta,u)=f(\theta,u)$

initialise $\theta$

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

sample $u\sim q(u|\theta)$

invoke involution $(\theta^{\prime},u^{\prime})=f(x,u)$

sample $\mathbf{y}\sim p(\mathbf{y}|\theta^{\prime})$

evaluate
$\alpha(\theta,\theta^{\prime})=\min\left\\{1,\frac{p(\theta^{\prime})\gamma(\mathbf{x}|\theta^{\prime})\gamma(\mathbf{y}|\theta)q(u^{\prime}|\theta^{\prime})}{p(\theta)\gamma(\mathbf{x}|\theta)\gamma(\mathbf{y}|\theta^{\prime})q(u|\theta)}\left|\frac{\partial
f(\theta,u)}{\partial(\theta,u)}\right|\right\\}$

$\theta_{i}=\begin{cases}\theta^{\prime}&\text{with probability
}\alpha(\theta,\theta^{\prime})\\\ \theta&\text{with probability
}1-\alpha(\theta,\theta^{\prime})\end{cases}$

$\theta\leftarrow\theta_{i}$

end for

Algorithm 1 Involutive exchange (iExchnage) algorithm

Let us first set the context. We have some data $\mathbf{x}$ which is assumed
to have been drawn via a model $p(\mathbf{x}|\theta)$, where $\theta$ denote
parameters, taking the following form

$\displaystyle
p(\mathbf{x}|\theta)=\frac{\gamma(\mathbf{x}|\theta)}{Z(\theta)}$ (18)

where $Z(\theta)=\int\gamma(\mathbf{x}|\theta)d\mathbf{x}$ denotes its
normalising constant, assumed to be intractable. If one is taking a Bayesian
approach to inference and has specified a prior $p(\theta)$, this leads to the
following posterior

$\displaystyle
p(\theta|\mathbf{x})=\frac{p(\mathbf{x}|\theta)p(\theta)}{p(\mathbf{x})}$ (19)

where $p(\mathbf{x})=\int p(\mathbf{x}|\theta)p(\theta)d\theta$ is the
marginal probability of the data, which in most cases is also intractable. Due
to these two elements of intractability, such posteriors are often referred to
as doubly-intractable (Murray et al., 2006). For example, the posteriors
resulting from both our SIS and SIM models are doubly-intractable.

A typical approach to circumvent the intractability present in Bayesian
posterior distributions is to use MCMC algorithms to sample from them, with
the Metropolis-Hastings (MH) algorithm being a prevalent choice. However, for
doubly-intractable posteriors, many standard MCMC algorithms are not feasible.
To illustrate, consider using the MH algorithm. Here, with $\theta$ the
current state and $q(\theta^{\prime}|\theta)$ some proposal density, in a
single iteration one would sample proposal $\theta^{\prime}$ from
$q(\theta^{\prime}|\theta)$ and accept this with the following probability

$\displaystyle\alpha(\theta,\theta^{\prime})$
$\displaystyle=\min\left\\{1,\frac{p(\theta^{\prime}|\mathbf{x})q(\theta|\theta^{\prime})}{p(\theta|\mathbf{x})q(\theta^{\prime}|\theta)}\right\\}$
(20)
$\displaystyle=\min\left\\{1,\frac{\gamma(\mathbf{x}|\theta^{\prime})/Z(\theta^{\prime})p(\theta^{\prime})q(\theta|\theta^{\prime})}{\gamma(\mathbf{x}|\theta)/Z(\theta)p(\theta)q(\theta^{\prime}|\theta)}\right\\},$

so that, starting from some initial state $\theta_{0}$ one obtains a sample
$\\{\theta_{i}\\}_{i=1}^{m}$ which is (approximately) distributed according to
$p(\theta|\mathbf{x})$. However, though the marginal probability
$p(\mathbf{x})$ cancels out in (20), the normalising constants $Z(\theta)$ and