# Nonparametric Two-Sample Test for Networks Using Joint Graphon Estimation

Benjamin Sischka and Göran Kauermann
Department of Statistics, Ludwig-Maximilians-Universität München  This work
was partially funded by the BERD@NFDI Consortium (https://www.berd-nfdi.de).

###### Abstract

This paper focuses on the comparison of networks on the basis of statistical
inference. For that purpose, we rely on smooth graphon models as a
nonparametric modeling strategy that is able to capture complex structural
patterns. The graphon itself can be viewed more broadly as density or
intensity function on networks, making the model a natural choice for
comparison purposes. Extending graphon estimation towards modeling multiple
networks simultaneously consequently provides substantial information about
the (dis-)similarity between networks. Fitting such a joint model—which can be
accomplished by applying an EM-type algorithm—provides a joint graphon
estimate plus a corresponding prediction of the node positions for each
network. In particular, it entails a generalized network alignment, where
nearby nodes play similar structural roles in their respective domains. Given
that, we construct a chi-squared test on equivalence of network structures.
Simulation studies and real-world examples support the applicability of our
network comparison strategy.

Keywords: Network Comparison, EM Algorithm, Gibbs Sampler, B-Spline
Regression, Chi-Squared Test

##  1 Introduction

The field of statistical modeling and analysis of complex networks has gained
strongly increasing interest over the last two decades. This is driven by the
fact that different types of systems can be reasonably formalized as
relationships between individuals or interactions between objects. Analyzing
such structures consequently allows to uncover and describe the phenomena that
affect these systems. Network-structured data arise in various fields, for
example, social and political sciences, economics, biology, neurosciences, and
many others. In this regard, a connectivity pattern between entities might
describe friendships among members of a social group (Eagle et al., 2009), the
trading between nations (Bhattacharya et al., 2008), interactions of proteins
(Schwikowski et al., 2000), or the functional coactivation within the human
brain (Bassett et al., 2018, Crossley et al., 2013).

In many situations, uncovering the underlying connectivity structure is not
the only concern but also the comparison of akin networks and the exploration
of potential differences. For example, this might be of interest in the
context of brain coactivation. Recently, a lot of work has been going on
investigating how the functional connectivity in the brain differs when people
are affected by cognitive disorders like Alzheimer’s disease or autism
spectrum disorder (Song et al., 2019, Subbaraju et al., 2017, Pascual-Belda et
al., 2018). Two such functional coactivation networks—resulting from
respectively averaging over the measurements of two different subject
groups—are illustrated in Figure 1.

Figure 1: Functional coactivation networks of the human brain. The illustrated
connectivity patterns result from averaging over multiple measurements for
subjects with autism spectrum disorder (left) and typical development (right).
Do these networks reveal a significant structural difference?

The posed research question in this context apparently also involves the
inquiry of whether a significant difference in the brain processes is
observable at all, which additionally might depend on the environmental
conditions like resting state, external stimuli, etc. More generally, this can
be phrased as a hypothesis test on structural equivalence of two networks.

To this end, we pursue constructing a model-based approach for network
comparison that allows for formal statistical testing. More precisely, we aim
to test whether two networks can be considered independent samples drawn from
the same probability distribution. This is apparently in itself a technically
difficult question since the two networks can have different sizes, and hence
the two distributions need to be somehow different. Therefore, it is crucial
that the applied distributional framework constitutes a rather universal
probability measure. In fact, this is not a trivial property, and many network
models entail conceptual issues that impede a direct comparison. For example,
in the Exponential Random Graph Model (Robins et al., 2007), a concrete model
parameterization has different implications for different network sizes,
making corresponding coefficient estimates hardly comparable. Hence, it is
necessary to rely on distributional models where the specification of size and
edge probabilities can be explicitly disentangled.

For such a comparative analysis of networks, we will demonstrate that Graphon
Models (Lovász and Szegedy, 2006, Diaconis and Janson, 2008) are a very useful
tool. First, the graphon model is very flexible and able to capture complex
network structures. Secondly, the graphon itself can be interpreted as
nonparametric density or intensity function on networks. Both together make
the graphon an overall characterizing network feature that, like nothing else,
uniquely covers the structure in a comprehensive way. Hence, the graphon
framework appears as a natural choice for comparison purposes. Lastly, the
model’s design fulfills the above requirement of decoupling the network’s
structure and size, allowing for modeling multiple networks simultaneously
(Navarro and Segarra, 2022).

The rest of the paper is organized as follows. In Section 2, we start with
reviewing methods from the network comparison literature. A formalization of
the test problem we want to tackle in this work is then concretely specified
and discussed in Section 3. The involved smooth graphon estimation in its
capacity as the joint modeling approach is formulated in Section 4. Based on
this joint graphon model, in Section 5, we develop a network comparison
strategy for testing on equivalence of the underlying structures. The general
applicability of the complete approach is demonstrated in Section 6, where we
consider its performance on simulated and real-world data. This involves the
method’s ability to uncover the underlying structure by joint graphon
estimation, as well as the qualification of the subsequent testing procedure.
The discussion and conclusion in Section 7 completes the paper.

##  2 Concepts for Network Comparison

When reviewing the literature on network comparison, it is worth noting that
all proposed strategies are naturally based on a specific concept of capturing
network structures. In general, the various approaches available for comparing
networks can be broadly distinguished according to whether they rely on a
descriptive or a model-based structural framework. Survey articles in this
field are given by Soundarajan et al. (2014), Yaveroğlu et al. (2015), Emmert-
Streib et al. (2016), and Tantardini et al. (2019). A more general perspective
on how complex data objects—such as adjacency matrices—might be compared is
pointed out by Marron and Alonso (2014). Consulting this compendium clearly
reveals a lack of model-based approaches in the context of network comparison.
This is specifically deficient since drawing statistical inference is only
possible under some kind of distributional assumption. The different
strategies for network comparison proposed in the literature—irrespective of
the capacity for drawing inference—are briefly reviewed hereafter.

Starting with approaches that are based on extracted network statistics, the
most intuitive strategy is probably to simply compare global characteristics
such as the clustering coefficient or the average path length (Newman, 2018,
pp. 364 ff., Butts, 2008, p. 31). However, this captures the overall network
structure only very poorly since completely differently structured networks
can apparently still possess the same global statistics. As a more advanced
approach, Wilson and Zhu (2008) consider the differences in the graph spectra,
see also Gera et al. (2018). Yet, for the spectrum, it is often unclear which
information in terms of local structural properties is extracted from the
network. In addition, spectral methods can be strongly affected by small
structural changes under specific circumstances. Taken together, such
approaches often ascribe too little importance to the attributes of interest,
leading to an over- or underrating of the structural dissimilarity at hand.

Another branch of the literature on descriptive network comparison relies on
the concept of graphlets, i.e. prespecified subgraph patterns that are assumed
to be sufficient for describing the present structure. Papers that, in one way
or another, consider differences in the frequencies of graphlets are, among
others, Pržulj et al. (2004), Pržulj (2007), Ali et al. (2014), and Faisal et
al. (2017). Since the counting procedure is rather complex for larger
graphlets, it is sort of a consensus to include only those that consist of no
more than five nodes. However, this seems somehow arbitrary and incomplete in
terms of capturing all structural aspects. Moreover, Yaveroğlu et al. (2014)
found high correlations among the graphlet-related statistics, including
complete redundancies. In contrast, model-based approaches specify the
complete distribution of frequencies over all kinds of subgraphs. A connection
between subgraph frequencies and the concrete specification of the graphon
model is exemplarily elaborated in Borgs et al. (2008), Bickel et al. (2011),
and Latouche and Robin (2016).

Overall, descriptive network statistics entail two general shortcomings.
First, it is very difficult to assess in which parts of the networks the key
differences are accommodated. To be precise, the nodes or edges (present or
absent) that contribute most to a quantified structural discrepancy can only
hardly be detected. Second, and more importantly, descriptive methods lack, by
nature, the ability to draw inference in the probabilistic sense.
Specifically, they provide no information on whether the found deviation
between networks is plausible to be ascribed to randomness or whether there is
a significant structural dissimilarity. Butts (2008) aims to overcome this
deficit by applying a simplistic conditional uniform graph distribution.

On the other hand, probabilistic models for network data allow to induce
distributions on network patterns that extend to desired distributional
assumptions on structural differences. As a consequence, these modeling
approaches might potentially serve as structural construct used for comparison
purposes. Yet, they possess individual conceptual shortcomings that often
impede a direct comparison. While, for example, the Latent Distance Model
(Hoff et al., 2002) does not provide any model-related key component to be
compared, coefficient estimates from the exponential random graph model are
not directly comparable across separated networks. The graphon model and the
Stochastic Blockmodel (Holland et al., 1983, Snijders and Nowicki, 1997 and
2001) suffer from identifiability issues (see e.g. Diaconis and Janson, 2008,
Thm. 7.1) that make a comparison of corresponding individual estimates
complicated. The latter model’s adaptivity is additionally highly dependent on
the choice of the number of blocks. Onnela et al. (2012) tackle this issue by
observing the networks’ complete disintegration processes, which they
subsequently summarize by the profiles of well-known network statistics.
Integrating over the profiles’ differences and applying principle component
analysis for summarization yields the final distance measure, which has been
demonstrated to provide reasonable results in terms of leading to a good
classification. However, to the best of our knowledge, there exists no (model-
based) nonparametric test on the equivalence of network structures.

In this paper, we aim to address this shortcoming, which we tackle by striking
new paths. As a general concept for this approach, we follow the intuition of
fitting a joint model to multiple networks simultaneously. For that purpose,
we resort to the smooth graphon model as an appropriate and very powerful
framework. Such a joint modeling strategy consequently circumvents the need
for post-hoc alignment of individual model fits and yields an outcome that
provides substantial information for comparison purposes. More precisely, it
allows for directly relating the networks at hand on the microscopic scale,
which, in the literature, is often referred to as “network alignment”
(Kuchaiev et al., 2010). However, as an essential distinction to classical
network alignment strategies, this method does not seek to find a node-wise
one-to-one mapping. Instead, it implies a mapping of local components, meaning
small fuzzy groups of nodes with similar structural roles in their respective
domains. Based on this network alignment, a structural comparison at the
microscopic level becomes possible. Aggregating local differences finally
enables to construct a test on structural equivalence of networks.

##  3 Notation and Formulation of the Test Problem

We consider the setting where two undirected networks of possibly different
sizes $N^{(1)}$ and $N^{(2)}$ have been observed. Let
$\boldsymbol{y}^{(g)}=[y_{ij}^{(g)}]_{i,j=1,\ldots,N^{(g)}}$ for $g=1,2$
denote the two respective adjacency matrices, where, for
$i,j=1,\ldots,N^{(g)}$, $y_{ij}^{(g)}=1$ if in network $g$ an edge between
nodes $i$ and $j$ exists and $y^{(g)}_{ij}=0$ otherwise. That specifically
means that $\boldsymbol{y}^{(g)}\in\\{0,1\\}^{N^{(g)}\times N^{(g)}}$. We
assume the networks to be undirected so that $y_{ij}^{(g)}=y_{ji}^{(g)}$.
Additionally, the diagonal elements are set to zero, i.e. $y_{ii}^{(g)}=0$,
reflecting the absence of self-loops. In general, we consider
$\boldsymbol{y}^{(g)}$ to be a realization of a random network
$\boldsymbol{Y}^{(g)}$ of size $N^{(g)}$ which is subject to probability mass
$\mathbb{P}(\boldsymbol{Y}^{(g)}=\boldsymbol{y}^{(g)}\,;N^{(g)})$. The
question we aim to tackle is whether $\boldsymbol{y}^{(1)}$ and
$\boldsymbol{y}^{(2)}$ are drawn from the same distribution. To suitably
specify such a distribution, we rely on the smooth graphon model. The data-
generating process is thereby as follows. Assume that we independently draw
uniformly distributed random variables

$U^{(g)}_{i}\sim\mbox{Uniform}[0,1]\quad\mbox{for }i=1,\ldots,N^{(g)}\mbox{
and }g=1,2.$ (1)

Conditional on $\boldsymbol{U}^{(g)}=(U^{(g)}_{1},\ldots,U_{N^{(g)}}^{(g)})$,
we then draw the edges i.i.d. through

$Y_{ij}^{(g)}\mid(\boldsymbol{U}^{(g)}=\boldsymbol{u}^{(g)})\sim\mbox{Binomial}(1,w^{(g)}(u^{(g)}_{i},u^{(g)}_{j}))\quad\mbox{for
}j>i$ (2)

with $\boldsymbol{u}^{(g)}=(u^{(g)}_{1},\ldots,\allowbreak
u_{N^{(g)}}^{(g)})\in[0,1]^{N^{(g)}}$ and under the setting of
$Y_{ij}^{(g)}\equiv Y_{ji}^{(g)}$ for $j<i$ and $Y_{ii}^{(g)}\equiv 0$. In
this modeling framework, the function $w^{(g)}:[0,1]^{2}\rightarrow[0,1]$,
which specifies the structural behavior of the emerging network, is called
graphon (see Lovász and Szegedy, 2006 and Diaconis and Janson, 2008). Here, in
particular, we assume $w^{(g)}(\cdot,\cdot)$ to be smooth according to some
Hölder or Lipschitz condition (cf. Wolfe and Olhede, 2013 or Chan and Airoldi,
2014). Relying on this data-generating process, the graphon-based probability
model can be defined through

$\boldsymbol{Y}^{(g)}\sim\mathbb{P}\big{(}\boldsymbol{Y}^{(g)}=\cdot\,;\,w^{(g)}(\cdot,\cdot),N^{(g)}\big{)}.$
(3)

In this distribution model, the network’s size and structure are apparently
dissociated, which therefore allows for a size-independent comparison of
underlying structures. Hence, our goal is to develop a statistical test on the
hypothesis

$H_{0}:\;w^{(1)}(\cdot,\cdot)\equiv w^{(2)}(\cdot,\cdot).$ (4)

In this context, we emphasize that data-generating process (2) is not unique
because it is invariant to permutations of $w^{(g)}(\cdot,\cdot)$, as
discussed in detail by Diaconis and Janson (2008, Sec. 7). Thus, the
formulation of $H_{0}$ needs to be understood from the perspective of
corresponding equivalence classes, implying that $w^{(1)}(\cdot,\cdot)$ and
$w^{(2)}(\cdot,\cdot)$ are rather viewed from a theoretical perspective.
Nonetheless, for the concrete implementation of the test procedure, we employ
a concrete representation of the two graphons. Specifically, under the
assumption of $H_{0}$ being true, we call the coinciding manifestation the
joint graphon. This can be formalized as

$w^{\text{joint}}(u,v):=w^{(1)}(u,v)=w^{(2)}(u,v)\quad\mbox{for all
}(u,v)^{\top}\in[0,1]^{2}.$

Since $w^{(1)}(\cdot,\cdot)$ and $w^{(2)}(\cdot,\cdot)$ are assumed to be
smooth—at least for one possible arrangement, and, in particular, the one we
consider here—, this also holds for $w^{\text{joint}}(\cdot,\cdot)$. Given
such a concrete representation of the joint graphon, the node position vectors
$\boldsymbol{u}^{(1)}$ and $\boldsymbol{u}^{(2)}$ referring to
$w^{\text{joint}}(\cdot,\cdot)$ then provide a specific type of network
alignment. This is what we utilize for a direct comparison of
$\boldsymbol{y}^{(1)}$ and $\boldsymbol{y}^{(2)}$. However, one typically
observes neither $\boldsymbol{u}^{(1)}$ and $\boldsymbol{u}^{(2)}$ nor
$w^{\text{joint}}(\cdot,\cdot)$. Thus, in order to achieve this alignment, we
first need to formulate an appropriate estimation procedure for the joint
graphon model.

##  4 EM-Based Joint Graphon Estimation

In this section, we present an iterative estimation procedure for the joint
smooth graphon model under the assumption that null hypothesis (4) is true. To
do so, we follow the EM-based estimation routine of Sischka and Kauermann
(2022a), extending it to the situation of two networks.

###  4.1 MCMC E-Step

Starting with the E-step of our iterative algorithm, we assume the joint
graphon $w^{\text{joint}}(\cdot,\cdot)$ to be known for the moment. Based on
that, the latent positions of the networks can be separately determined using
MCMC techniques. To be precise, we apply Gibbs sampling by formulating the
full conditional distribution of $U_{i}^{(g)}$ through

$f(u_{i}^{(g)}\mid
u_{1}^{(g)},\ldots,u_{i-1}^{(g)},u_{i+1}^{(g)},\ldots,u_{N^{(g)}}^{(g)},\boldsymbol{y}^{(g)})\\\
\propto\prod_{j\neq
i}w^{\text{joint}}(u_{i}^{(g)},u_{j}^{(g)})^{y_{ij}^{(g)}}[1-w^{\text{joint}}(u_{i}^{(g)},u_{j}^{(g)})]^{1-{y_{ij}^{(g)}}}$
(5)

for all $i=1,\ldots,N^{(g)}$ and $g=1,2$. Details on the concrete
implementation of the Gibbs sampler are given in Section A of the Appendix.
The resulting MCMC sequence (after cutting the burn-in period and appropriate
thinning) then reflects the joint conditional distribution
$f(\boldsymbol{u}^{(g)}\mid\boldsymbol{y}^{(g)})$. Thus, the marginal
conditional means of the latent positions, i.e.
$\mathbb{E}(U_{i}^{(g)}\mid\boldsymbol{Y}^{(g)}=\boldsymbol{y}^{(g)})$ for
$i=1,\ldots,N^{(g)}$, can be approximated by taking the mean over the MCMC
samples, which we denote by
$\bar{\boldsymbol{u}}^{(g)}=(\bar{u}_{1}^{(g)},\ldots,\bar{u}_{N^{(g)}}^{(g)})$.
This posterior mean vector, however, requires further adjustments due to
additional identifiability issues which cannot be coped with the standard EM-
type algorithm. To illustrate this, let model assumption (1) be more relaxed
in the sense that the $U_{i}^{(g)}$’s might follow any continuous distribution
$F^{(g)}(\cdot)$. Under this configuration, the model ($F^{(g)}(\cdot)$,
$w^{(g)}(\cdot,\cdot)$) is equivalent to any other model
(${F^{(g)}}^{\prime}(\cdot)$, ${w^{(g)}}^{\prime}(\cdot,\cdot)$) constructed
through

${F^{(g)}}^{\prime}(u^{\prime}):=F^{(g)}(\varphi(u^{\prime}))\quad\mbox{and}\quad{w^{(g)}}^{\prime}(u^{\prime},v^{\prime}):=w^{(g)}(\varphi(u^{\prime}),\varphi(v^{\prime}))$

with $\varphi:[0,1]\rightarrow[0,1]$ being a strictly increasing continuous
function (that is, in contrast to Diaconis and Janson, 2008, Sec. 7, not
measure-preserving). Specifically, that means

$\mathbb{P}\big{(}\boldsymbol{Y}^{(g)}=\cdot\,;\,F^{(g)}(\cdot),\allowbreak
w^{(g)}(\cdot,\cdot),\allowbreak
N^{(g)}\big{)}\equiv\mathbb{P}\big{(}\boldsymbol{Y}^{(g)}=\cdot\,;\,{F^{(g)}}^{\prime}(\cdot),\allowbreak{w^{(g)}}^{\prime}(\cdot,\cdot),\allowbreak
N^{(g)}\big{)}$

for all $N^{(g)}\geq 2$. As a matter of conception, this issue cannot be
solved by the EM algorithm since it aims at specifying a model that adapts
optimally to the given data instead of perfectly recovering the underlying
model structure. Consequently, the EM approach is not able to distinguish
between the two conceptually equivalent model specifications
($F^{(g)}(\cdot)$, $w^{(g)}(\cdot,\cdot)$) and (${F^{(g)}}^{\prime}(\cdot)$,
${w^{(g)}}^{\prime}(\cdot,\cdot)$). Nonetheless, this identifiability issue
can simply be tackled by adjusting $\bar{\boldsymbol{u}}^{(g)}$ before
estimating the graphon in the M-step. To do so, we just impose that the
inferred node positions follow an ideal sample drawn from the standard uniform
distribution. That is, we set

$\hat{u}_{i}^{(g)}=\frac{\operatorname{rank}(\bar{u}_{i}^{(g)})}{N^{(g)}+1},$

where $\operatorname{rank}(\bar{u}_{i}^{(g)})$ is the rank from smallest to
largest of element $\bar{u}_{i}^{(g)}$ within $\bar{\boldsymbol{u}}^{(g)}$. In
this context, note that the values $i/(N^{(g)}+1)$ with $i=1,\ldots,N^{(g)}$
represent the expectations of $N^{(g)}$ ordered random variables that are
independently drawn from the standard uniform distribution. As a result, with
$\hat{\boldsymbol{u}}^{(g)}=(\hat{u}_{1}^{(g)},\ldots,\hat{u}_{N^{(g)}}^{(g)})$
we obtain a plausible realization of the node positions of network $g$.
Apparently, this relies on the current joint graphon estimate
$\hat{w}^{\text{joint}}(\cdot,\cdot)$, which is applied as substitute in
conditional distribution (5). In the next step, we formulate the procedure for
updating $\hat{w}^{\text{joint}}(\cdot,\cdot)$ given
$\hat{\boldsymbol{u}}^{(1)}$ and $\hat{\boldsymbol{u}}^{(2)}$.

###  4.2 Spline-Based M-Step

For a semiparametric estimation of the joint smooth graphon, we choose a
linear B-spline regression approach. To this end, we assume the joint graphon
to be approximated through

$w_{\boldsymbol{\theta}}^{\text{joint}}(u,v)=\boldsymbol{B}(u,v)\,\boldsymbol{\theta}=\left[\boldsymbol{B}(u)\otimes\boldsymbol{B}(v)\right]\boldsymbol{\theta},$

where $\otimes$ is the Kronecker product,
$\boldsymbol{B}(u)\in\mathbb{R}^{1\times L}$ is a linear B-spline basis on
$[0,1]$, normalized to have a maximum value of one, and
$\boldsymbol{\theta}\in\mathbb{R}^{L^{2}}$ is the parameter vector to be
estimated. The inner B-spline knots are specified as lying equidistantly on a
regular 2D grid within $[0,1]^{2}$, where $\boldsymbol{\theta}$ is indexed
accordingly through
$\boldsymbol{\theta}=\left(\theta_{11},\ldots,\theta_{1L},\allowbreak\theta_{21},\ldots,\theta_{LL}\right)^{\top}$.
Based on this representation and given the node positions
$\hat{\boldsymbol{u}}^{(1)}$ and $\hat{\boldsymbol{u}}^{(2)}$, we formulate
the marginal log-likelihood over both networks as

$\ell(\boldsymbol{\theta})=\sum_{g}\sum\limits_{\begin{subarray}{c}i,j\\\
j\neq
i\end{subarray}}\left[y_{ij}^{(g)}\,\log\left(\boldsymbol{B}_{ij}^{(g)}\boldsymbol{\theta}\right)+\left(1-y_{ij}^{(g)}\right)\,\log\left(1-\boldsymbol{B}_{ij}^{(g)}\boldsymbol{\theta}\right)\right],$
(6)

where
$\boldsymbol{B}_{ij}^{(g)}=\boldsymbol{B}(\hat{u}_{i}^{(g)})\otimes\boldsymbol{B}(\hat{u}_{j}^{(g)})$.
Furthermore, through standard calculations, we are able to derive the score
function $\boldsymbol{s}(\boldsymbol{\theta})$ and the Fisher information
$\boldsymbol{F}(\boldsymbol{\theta})$, as demonstrated in Section B of the
Appendix. Fisher scoring can then be used to maximize
$\ell(\boldsymbol{\theta})$. In addition, we include side constraints to
ensure that $w_{\boldsymbol{\theta}}^{\text{joint}}(\cdot,\cdot)$ is bounded
to $[0,1]$ and symmetric. In the linear B-spline setting, this means
restricting the parameters by the conditions

$\theta_{kl}\geq 0\;,\quad\theta_{kl}\leq
1\;,\quad\mbox{and}\quad\theta_{kl}-\theta_{lk}=0$

for all $l>k$. Apparently, all three conditions are of linear form and thus
can be written in matrix format. Taken together, the Fisher scoring becomes a
quadratic programming problem that can be solved using standard software (see
Andersen et al., 2016 or Turlach and Weingessel, 2013).

Moreover, we intend to add penalization on the B-spline estimate. As outlined
in Eilers and Marx (1996) and Ruppert et al. (2003, 2009), penalized spline
estimation under the setting of a rather large spline basis yields a
preferable outcome since it guarantees a functional fit that covers the data
adequately but is still smooth. Thus, this approach enables to precisely
capture the underlying structure while avoiding overfitting. To realize this,
we add a first-order penalty, meaning that “neighboring” elements of
$\boldsymbol{\theta}$ get penalized. For the log-likelihood, the score
function, and the Fisher information, this leads to the penalized versions in
the form of

$\begin{gathered}\ell_{\text{p}}(\boldsymbol{\theta},\lambda)=\ell(\boldsymbol{\theta})-\frac{1}{2}\lambda\boldsymbol{\theta}^{\top}\boldsymbol{P}\boldsymbol{\theta}\,,\quad\boldsymbol{s}_{\text{p}}(\boldsymbol{\theta},\lambda)=\boldsymbol{s}(\boldsymbol{\theta})-\lambda\boldsymbol{P}\boldsymbol{\theta}\,,\\\
\text{and}\quad\boldsymbol{F}_{\text{p}}(\boldsymbol{\theta},\lambda)=\boldsymbol{F}(\boldsymbol{\theta})+\lambda\boldsymbol{P},\end{gathered}$
(7)

respectively, where $\boldsymbol{P}$ is a penalization matrix of appropriate
shape (see Section B of the Appendix). For an adequate choice of the penalty
parameter $\lambda$ in the two/dimensional spline regression, we follow
Kauermann et al. (2013) and apply the corrected Akaike Information Criterion
($\operatorname{AIC}_{\text{c}}$, see Hurvich and Tsai, 1989 and Burnham and
Anderson, 2002). This is defined as

$\operatorname{AIC}_{\text{c}}(\lambda)=-2\,\ell(\hat{\boldsymbol{\theta}}_{\text{p}})+2\,\operatorname{df}(\lambda)+\frac{2\,\operatorname{df}(\lambda)[\operatorname{df}(\lambda)+1]}{N(N-1)-\operatorname{df}(\lambda)-1},$

where $\hat{\boldsymbol{\theta}}_{\text{p}}$ is the corresponding penalized
parameter estimate and $\operatorname{df}(\lambda)$ specifies the degrees of
freedom of the penalized B-spline function. More precisely, according to Wood
(2017, pp. 211 ff.), the latter is defined trough

$\operatorname{df}(\lambda)=\operatorname{tr}\left\\{\boldsymbol{F}_{\text{p}}^{-1}(\hat{\boldsymbol{\theta}}_{\text{p}},\lambda)\,\boldsymbol{F}(\hat{\boldsymbol{\theta}}_{\text{p}})\right\\}$

with $\operatorname{tr}\\{\cdot\\}$ being the trace of a matrix. A numerical
optimization of the corrected $\operatorname{AIC}$ with respect to $\lambda$
concludes the estimation of $\boldsymbol{\theta}$, resulting in the eventual
estimate $\hat{w}^{\text{joint}}(\cdot,\cdot)$ of the current M-step.

Finally, the EM-type estimation procedure described above—meaning the
consecutive repetition of the E- and M-step until convergence is
achieved—allows us to adequately estimate both the joint graphon
$w^{\text{joint}}(\cdot,\cdot)$ and the corresponding node positions
$\boldsymbol{u}^{(1)}$ and $\boldsymbol{u}^{(2)}$ of the two networks. Based
on these results, we are now able to formulate an appropriate test procedure.

##  5 Two-Sample Test on Network Structures

Returning to the test problem raised in Section 3, we now develop a
statistical test procedure on hypothesis (4), i.e. whether
$\boldsymbol{y}^{(1)}$ and $\boldsymbol{y}^{(2)}$ are drawn from the same
distribution. To do so, we utilize the network alignment resulting from the
(inferred) joint smooth graphon model. More precisely, we exploit the fact
that two edge variables $Y_{i_{1}j_{1}}^{(1)}$ and $Y_{i_{2}j_{2}}^{(2)}$ that
have nearby positions—i.e. for which the distance between
$(u_{i_{1}}^{(1)},u_{j_{1}}^{(1)})^{\top}$ and
$(u_{i_{2}}^{(2)},u_{j_{2}}^{(2)})^{\top}$ is small—possess similar
probabilities to form a connection. In a more formalized way, this means that,
from
$\|(u_{i_{1}}^{(1)},u_{j_{1}}^{(1)})^{\top}-(u_{i_{2}}^{(2)},u_{j_{2}}^{(2)})^{\top}\|\approx
0$, it follows that

$\mathbb{P}(Y_{i_{1}j_{1}}^{(1)}=1\mid
U_{i_{1}}^{(1)}=u_{i_{1}}^{(1)},U_{j_{1}}^{(1)}=u_{j_{1}}^{(1)})\\\
\approx\mathbb{P}(Y_{i_{2}j_{2}}^{(2)}=1\mid
U_{i_{2}}^{(2)}=u_{i_{2}}^{(2)},U_{j_{2}}^{(2)}=u_{j_{2}}^{(2)}),$

where $\|\cdot\|$ is the Euclidean distance. Following this intuition, we
divide the unit square into small segments and compare between networks the
ratio of present versus absent edges occurring in these segments (see Figure 2
for an exemplary division).

Figure 2: Dividing the unit square as domain of the graphon model into small
segments for comparing network structure on a microscopic level. Left:
division of $w^{\text{joint}}(\cdot,\cdot)$ into approximately piecewise-
constant rectangles. Middle and right: edge positions
$(u_{i}^{(g)},u_{j}^{(g)})^{\top}$ of two simulated networks with respect to
$w^{\text{joint}}(\cdot,\cdot)$; weakly colored crosses and intensively
colored circles represent absent and present edges, respectively. The two
networks can be compared by pairwise contrasting the edge proportions within
the labeled rectangles.

For that purpose, we choose a suitable $K\in\mathbb{N}$, specify a
corresponding boundary sequence $a_{0}=0<a_{1}<\ldots<a_{K}=1$, and define the
following two quantities for $l,k=1,\ldots,K$ with $l\geq k$:

$\displaystyle\begin{split}d_{kl}^{(g)}&=\sum\limits_{\begin{subarray}{c}i,j\\\
j>i\end{subarray}}y_{ij}^{(g)}\mathds{1}_{\\{u_{i}^{(g)}\in[a_{k-1},a_{k})\\}}\mathds{1}_{\\{u_{j}^{(g)}\in[a_{l-1},a_{l})\\}}\\\
m_{kl}^{(g)}&=\sum\limits_{\begin{subarray}{c}i,j\\\
j>i\end{subarray}}\mathds{1}_{\\{u_{i}^{(g)}\in[a_{k-1},a_{k})\\}}\mathds{1}_{\\{u_{j}^{(g)}\in[a_{l-1},a_{l})\\}}.\end{split}$
(8)

This means $d_{kl}^{(g)}$ and $m_{kl}^{(g)}$ represent the number of present
($y_{ij}^{(g)}=1$) and a priori potential ($y_{ij}^{(g)}\in\\{0,1\\}$) edges
of network $g$, respectively, within the constructed rectangle
$[a_{k-1},a_{k})\times[a_{l-1},a_{l})$. The corresponding cross-network counts
can be calculated by $d_{kl}=d_{kl}^{(1)}+d_{kl}^{(2)}$ and
$m_{kl}=m_{kl}^{(1)}+m_{kl}^{(2)}$. Since $w^{\text{joint}}(\cdot,\cdot)$ is
smooth, we further assume that the induced probability on edge variables
within $[a_{k-1},a_{k})\times[a_{l-1},a_{l})$ is approximately constant. That
allows for putting the observed ratios between present and absent edges in
direct relation. In this light, we formulate the following contingency table
to keep track of homogeneity between the networks within rectangle $(k,l)$:

$d_{kl}^{(1)}$ | $d_{kl}^{(2)}$ | $d_{kl}$
---|---|---
$m_{kl}^{(1)}-d_{kl}^{(1)}$ | $m_{kl}^{(2)}-d_{kl}^{(2)}$ | $m_{kl}-d_{kl}$
$m_{kl}^{(1)}$ | $m_{kl}^{(2)}$ | $m_{kl}$

Apparently, if $H_{0}$ is assumed to be true, we would expect the proportions
of present edges, $d_{kl}^{(1)}/m_{kl}^{(1)}$ and $d_{kl}^{(2)}/m_{kl}^{(2)}$,
to be similar. This can be assessed by contrasting the observed numbers of
edges with their expectations conditional on the given margin totals, which is
in line with the construction of Fisher’s exact test on $2\times 2$
contingency tables. In this regard, the theoretical random counterpart of
$d_{kl}^{(1)}$ can be defined as

$\displaystyle D_{kl}^{(1)}$
$\displaystyle=\sum\limits_{\begin{subarray}{c}i,j\\\
j>i\end{subarray}}Y_{ij}^{(1)}\mathds{1}_{\\{u_{i}^{(1)}\in[a_{k-1},a_{k})\\}}\mathds{1}_{\\{u_{j}^{(1)}\in[a_{l-1},a_{l})\\}},$

for which under $H_{0}$ it approximately holds that

$\begin{gathered}D_{kl}^{(1)}\mid
d_{kl}\sim\mbox{Hyp}\left(m_{kl},d_{kl},m_{kl}^{(1)}\right)\quad\mbox{with}\quad
E_{kl}^{(1)}:=\mathbb{E}(D_{kl}^{(1)}\mid
d_{kl})=m_{kl}^{(1)}\frac{d_{kl}}{m_{kl}}\\\ \mbox{and}\quad
V_{kl}^{(1)}:=\mathbb{V}(D_{kl}^{(1)}\mid
d_{kl})=m_{kl}^{(1)}\frac{d_{kl}}{m_{kl}}\frac{m_{kl}-d_{kl}}{m_{kl}}\frac{m_{kl}-m_{kl}^{(1)}}{m_{kl}-1}.\end{gathered}$
(9)

Based on these specifications, we define our test statistic as

$T=\sum\limits_{\begin{subarray}{c}k,l\\\ l\geq
k\end{subarray}}\frac{\left(D_{kl}^{(1)}-E_{kl}^{(1)}\right)^{2}}{V_{kl}^{(1)}}\quad\mbox{with
realization}\quad t=\sum\limits_{\begin{subarray}{c}k,l\\\ l\geq
k\end{subarray}}\frac{\left(d_{kl}^{(1)}-E_{kl}^{(1)}\right)^{2}}{V_{kl}^{(1)}}.$
(10)

Note that we only include the quantities of the first network due to the
symmetry of the hypergeometric distribution. Moreover, summands for which
$V_{kl}^{(1)}=0$—resulting from $m_{kl}^{(1)}$, $m_{kl}-m_{kl}^{(1)}$,
$d_{kl}$, or $m_{kl}-d_{kl}$ being zero—carry no relevant information and thus
can simply be omitted from the calculation. In contrast, if $m_{kl}^{(1)}$ is
large, $m_{kl}$ and $d_{kl}$ are large compared to $m_{kl}^{(1)}$, and
$d_{kl}/m_{kl}$ is not close to zero or one, then $D_{kl}^{(1)}$ is known to
be approximately normally distributed. Given that, we can conclude that

$T\stackrel{{\scriptstyle\text{a}}}{{\sim}}\chi^{2}_{K(K+1)/2}$ (11)

since, in this scenario, the test statistic is essentially the sum of squared
(conditionally) independent random variables that approximately follow a
standard normal distribution. If the latter condition does not apply, and
assumption (11) is not reasonable to hold, we still can simulate a sample of
the theoretical distribution by drawing $D_{kl}^{(1)}\mid d_{kl}$ according to
(9) and calculating $T$ as in (10). In both cases, we can easily derive a
critical value $c_{1-\alpha}$ to be compared with the realization $t$ of the
test statistic. To do so, we pick the corresponding $(1-\alpha)$/quantile of
either the theoretical distribution $\chi^{2}_{K(K+1)/2}$ or the simulated
sample. Finally, we reject null hypothesis (4) at the significance level of
$\alpha$ if $t>c_{1-\alpha}$. The choice of an appropriate $K$ applied for
these calculations is discussed in Section C of the Appendix. Note that
altogether the presented test procedure follows a conception similar to the
one underlying the log-rank test for time-to-event data.

Apparently, when conducting the test procedure on real-world networks, we
obtain the joint graphon and the corresponding alignment of the networks by
applying the estimation procedure described in Section 4. In the end, this
enables us to appropriately approximate test statistic (10). In this context,
it is important to consider the general behavior of the joint graphon
estimation under the alternative, that is, if hypothesis (4) does not hold. We
stress that the intuition of the estimation procedure is to align the two
networks as well as possible with respect to some suitable joint graphon
model. Consequently, the expectation of $T$ will be higher the more the true
graphons $w^{(1)}(\cdot,\cdot)$ and $w^{(2)}(\cdot,\cdot)$ differ after
“optimal” alignment. This clearly implies that the power of our test is higher
for instances that deviate more strongly from the null hypothesis.

##  6 Applications

In this section, we showcase the applicability of the joint graphon estimation
routine and the subsequent testing procedure. To give a comprehensive insight,
this comprises both simulated and real-world networks. For an optimal
estimation result and to best approximate test statistic (10), we repeat the
estimation and testing procedure several times. In a modeling-oriented
context, we would then typically pick the outcome with the lowest corrected
$\operatorname{AIC}$. However, since here the focus is on the statistical
testing aspect, we choose the estimation result which leads to the highest
$p$-value, assuming that this provides an optimal lower bound for the outcome
under the (potentially existing) oracle network alignment.

###  6.1 Simulation Studies

####  6.1.1 Exemplary Application to Synthetic Data

To demonstrate the general capability of the joint graphon estimation and the
performance of the subsequent testing procedure, we consider the graphon in
the top left plot of Figure 3.

Figure 3: Joint graphon estimation for simulated networks with subsequent
testing on equivalence of the underlying distribution models. The top row
shows the true and the jointly estimated graphon on the left and right,
respectively. The realizations of the terms of test statistic (10),
representing the dissimilarities of the two networks per rectangle, are
visualized at the bottom left, where $m_{kl}^{(g)}\geq 100$ for $k\neq l$ and
$\geq 45$ otherwise. The final result of the test statistic (black solid
vertical line) as well as its distribution under $H_{0}$ are illustrated at
the bottom right, where the black solid step function and the blue dashed
curve depict the simulated and the asymptotic chi-squared distribution,
respectively. The red dashed vertical lines visualize the critical values at a
significance level of $5\%$, derived from the simulated (upper line) and the
asymptotic distribution (lower line).

Its formation is inspired by and can be interpreted as a stochastic blockmodel
with smooth transitions between communities. Based on this ground-truth model
specification, we simulate two networks with $N^{(1)}=200$ and $N^{(2)}=300$
by making use of data-generating process (2). To recover the underlying
structure, we then apply the presented EM-type algorithm, where, for
initialization, we make use of an uninformative random node positioning. After
several iterations, we achieve the reasonable joint graphon estimate at the
top right, which fully captures the structure of the ground-truth graphon.
Relying on the accompanying estimated node positions, we subsequently conduct
the testing procedure on whether the underlying distributions are equivalent.
To this end, we start with calculating the rectangle-wise differences
according to the construction of test statistic (10). The results are depicted
as a heat map at the bottom left plot of Figure 3. This reveals that the
difference in the local edge density is rather low to moderate in most
rectangles, whereas it is distinctly higher in a few others. Aggregating these
differences yields a test statistic of $203.2$ as depicted by the black solid
vertical line at the bottom right. Contrasting this result with the simulated
$95\%$ quantile of the distribution of $T$ under $H_{0}$ as the critical value
(red dashed vertical line) yields no rejection. Hence, the underlying
distributions of the two networks do not significantly differ with respect to
a significance level of $5\%$. As a final remark with regard to the bottom
right plot, the simulated distribution of $T$ (black solid step function) and
its theoretical approximation (blue dashed curve)—both relying on the
assumption of $H_{0}$ being true—are very close to one another. Consequently,
they also provide very similar critical values, namely $243.6$ and $244.8$,
respectively. This demonstrates that asymptotic distribution (11) represents a
good approximation.

####  6.1.2 Performance Analysis under $H_{0}$

To evaluate the performance of the testing procedure in this example more
profoundly, we repeat the above proceeding $400$ times, with newly simulated
networks in each trial (remaining with $N^{(1)}=200$ and $N^{(2)}=300$). Note
that we run the estimation procedure always ten times (with different random
node positions as varying initialization) and finally pick the highest
$p$-value as the actual result for the given network pair. These repetitions
already provide a broad insight into the method’s performance under the given
setting. An even more extensive evaluation becomes possible when, in contrast
to the proceeding above, the testing procedure is performed on the basis of
the oracle node positions. This allows us to dramatically reduce the
computational burden since it releases us from the preceding (computationally
expensive) model estimation. As a consequence, we are able to increase the
number of conducted tests to $10,000$. From these two repetition studies
(using either $\hat{\boldsymbol{u}}^{(g)}$ or $\boldsymbol{u}^{(g)}$), we
obtain rejection rates of $6.5\%$ and $6.15\%$ under the estimated and oracle
node positioning, respectively. That means the test is slightly overconfident
relative to the nominal significance level of $5\%$. The top row of Figure 4

Figure 4: Performance of the testing procedure with regard to the resulting
$p$-value; results are simulation-based. Top: empirical distribution of the
$p$-value under $H_{0}$, illustrated as density and cumulative distribution
function on the left and right, respectively. The black dashed lines
illustrate the desired distributional behavior of an optimal test. Number of
repetitions for estimated / oracle node positions: $400$ / $10,000$. Bottom:
distribution of the $p$-value under $H_{1}$ and the usage of oracle node
positions (in box plot format); based on $1,000$ repetitions each. The x-axis
illustrates different settings according to formulation (12) (higher value of
$\gamma$ implies stronger deviation from $H_{0}$). The black dashed horizontal
line represents the $5\%$ significance level, and the orange curve illustrates
the corresponding power.

shows additionally the empirical distributions of the observed $p$-values,
illustrated as densities (left) and cumulative distribution functions (right).
In accordance with the mildly inflated rejection rates, this exhibits a slight
tendency to underestimate the $p$-value, i.e. interpreting the discrepancy as
too high in distributional terms.

####  6.1.3 Performance Analysis under $H_{1}$

Conclusively, we are interested in evaluating the test performance under a
false null hypothesis, which apparently requires formulating a suitable
alternative. To this end, we “shrink” the heterogeneity within the graphon
such that the present structure becomes less pronounced. The resulting graphon
specification consequently tends more towards an Erdős–Rényi model, with the
global density remaining unchanged. To be precise, based on
$w^{(1)}(\cdot,\cdot)$, we formulate

$\displaystyle w^{(2)}(u,v):=(1-\gamma)\,w^{(1)}(u,v)+\gamma\,\bar{w}^{(1)}$
(12)

with $\gamma\in[0,1]$ and $\bar{w}^{(1)}=\iint
w^{(1)}(u,v)\mathop{}\\!\mathrm{d}u\mathop{}\\!\mathrm{d}v$. Apparently,
increasing the mixing parameter $\gamma$ leads to a stronger deviation from
$H_{0}$. At the same time, this setting guarantees an optimal alignment of
$w^{(1)}(\cdot,\cdot)$ and $w^{(2)}(\cdot,\cdot)$, meaning that there exists
no rearrangement of $w^{(2)}(\cdot,\cdot)$ that is closer to
$w^{(1)}(\cdot,\cdot)$ than specification (12). For this experiment, we again
choose $N^{(1)}=200$ and $N^{(2)}=300$. Moreover, here we rely exclusively on
the oracle node positions. This provides a lower bound of the power since the
rejection rate can be expected to be higher when using estimated node
positions instead (cf. previous analysis under $H_{0}$). The results for this
setup are presented in the bottom plot of Figure 4, where the distribution of
the resulting $p$-value is illustrated for different settings of $\gamma$. The
orange curve additionally visualizes the resulting power, i.e. the proportion
of cases with $p<0.05$. These results clearly show that the probability of
detecting the false null hypothesis monotonically increases as the parameter
$\gamma$ gets larger, which underpins the appropriateness of our test
procedure.

Overall, the obtained simulation results demonstrate that the elaborated
estimation and testing procedure yields reasonable results for assessing
structural differences between networks. Building upon these findings, we next
want to investigate the method’s performance on real-world data.

###  6.2 Real-World Examples

####  6.2.1 Facebook Ego Networks

As a first real-world example, we consider two Facebook ego networks which
have been assembled by Leskovec and McAuley (2012) and are publicly available
on the Stanford Large Network Dataset Collection (Leskovec and Krevl, 2014).
The two ego networks consist of $333$ and $168$ individuals, respectively,
where the ego nodes are not included. An illustration of these networks is
given in the top row of Figure 5,

Figure 5: Comparison of two Facebook ego networks. Top: illustration of
networks with coloring referring to estimated node positions. Middle: ordered
adjacency matrices divided into blockwise segments. Bottom left: segment-wise
differences between the two networks with $m_{kl}^{(g)}\geq 100$ for $k\neq l$
and $\geq 45$ otherwise; gray rectangles do not contain any observed edges
($d_{kl}=0$) and thus provide no information. Bottom right: realization of
test statistic (black solid vertical line) plus corresponding distribution
under $H_{0}$ (black solid step function and blue dashed curve represent
simulated and asymptotic chi-squared distribution, respectively); critical
values derived from the two types of distributions are represented by the
upper and the lower red dashed vertical line.

with nodes being colored according to resulting positions. In both networks,
these final estimated node positions appear to be in line with the given
network structure in terms of reflecting the nodes’ embedding within the
network. Moreover, they seem to be aligned across networks. For example, in
both networks, the reddish nodes represent the rather centric individuals,
whereas the nodes from the dark blue spectrum constitute a moderately
interconnected branch that is more detached from the rest of the network.
However, the segment-wise differences depicted at the bottom left in Figure 5
exhibit some severe deviations. This can be clearly traced back to the
blockwise division of the adjacency matrices as it results from partitioning
the domain of edge positions (middle row). The aggregated differences
ultimately result in a test statistic that is far from the sector of plausible
values under the null hypothesis, as illustrated at the bottom right.
Consequently, we can conclude that the structural behavior in the two networks
differs significantly.

####  6.2.2 Human Brain Functional Coactivation Networks

In the second real-world application, we are concerned with differences in the
human brain coactivation structure. To be precise, we compare two types of
individuals, one with autism spectrum disorder (ASD) and the other with
typical development (TD). In particular, we are interested in whether the
functional connectivity within the brain significantly differs between these
two groups (cf. the introductory example from Figure 1). For this analysis, we
use resting-state functional magnetic resonance imaging data from the Autism
Brain Imaging Data Exchange project (ABIDE I, 2013). More specifically, we
employ preprocessed data provided by the Preprocessed Connectomes Project
platform (PCP, 2015). Based on these person-specific datasets, we first
calculate correlations between brain regions with respect to concurrent
activation over time. Aggregating the results of participants from the same
clinical group and employing an appropriate threshold finally yields the
network-structured data which we aim to compare. To be precise, by performing
the described preprocessing, we achieve for both groups, ASD and TD, a global
connectivity pattern between $116$ prespecified relevant brain regions. Note
that these regions are the same for both groups, which is why this could also
be viewed as a comparison task under known node correspondence. However, we
emphasize that in neurosciences, the transfer of competencies between brain
regions is a well-known phenomenon, wherefore the general functional
connectivity structure might be of greater relevance than the functional
connection between specific regions. Further details on the acquisition and
adequate transformation of the data are provided in Section D of the Appendix.

For analyzing the differences in the brain coactivation structure between the
two diagnostic groups, we again start with appropriately aligning the two
networks. This is apparently done by employing the joint graphon estimation
routine. The resulting node positions in relation to the embedding of nodes
within the networks are illustrated in the top row of Figure 6.

Figure 6: Comparison of functional coactivation in the human brain between
groups of subjects with autism spectrum disorder and with typical development.
The top row shows the networks of the ASD and the TD group on the left- and
right-hand side, respectively. All illustration aspects are equivalent to the
representation in Figure 5. The number of nodes per rectangle is again given
by $m_{kl}^{(g)}\geq 100$ for $k\neq l$ and $\geq 45$ otherwise, where N/A’s
in the blockwise differences result from $d_{kl}$ or $m_{kl}-d_{kl}$ being
zero.

Again this reveals a plausible allocation of the nodes in the joint graphon
model. The structural evolvement can further be evaluated by consulting the
correspondingly ordered adjacency matrices (see middle row), where the red
dashed lines represent the blockwise division resulting from the assignment of
the edge positions to the rectangles in $[0,1]$. At first view, the formed
structure looks quite similar in both matrices. Yet, on closer inspection,
some blocks can be found where the density seems considerably different. This
is also observed in the rectangle-wise differences depicted in the bottom left
plot. In the end, aggregating these differences leads to rejecting the null
hypothesis, as represented at the bottom right. To be precise, this test
decision is based on a resulting $p$-value of $0.013$ (with reference to the
simulated distribution under $H_{0}$).

Given that this outcome does not support an utterly unambiguous decision of
the conducted test procedure, one might additionally be interested in the
nature of the inferred differences. To address this, in Section
LABEL:sec:brainLocDiff of the Supplementary Material, we localize different
behavior between the networks on the microscopic scale. Note that this can
also be derived more or less directly from the joint graphon model. Moreover,
for comparison reasons, we repeated the above analysis for two randomly
selected disjoint subgroups of the TD group. According to the results
illustrated in Figure LABEL:fig:realWorld2c of the Supplementary Material, in
this scenario, we do not observe a significant overall deviation. This further
underlines the findings about the dissimilarity between the ASD and the TD
group.

##  7 Discussion and Conclusion

In the network comparison literature, the task of drawing statistical
inference appears to be an open challenge up to now. We addressed this
shortcoming in this paper by developing a nonparametric test on the
equivalence of network structures. To do so, we utilized the smooth graphon
model as a powerful tool for both describing and modeling the structure in
complex networks. More precisely, extending previous estimation approaches
towards a joint modeling framework allowed us to formulate a more generalized
network alignment. Given that, local structure comparison can be carried out
to uncover differences on the microscopic scale. Adequately aggregating these
local differences finally enables to construct an appropriate nonparametric
testing procedure on network data. Applying this comparison strategy to
simulated and real-world networks clearly demonstrated its general
applicability.

As outlined before, a crucial point for the proposed approach to work is the
graphon model’s property of decoupling structure and size. In the same line,
one could think of further decoupling the global density by following the
approach of Bickel and Chen (2009), i.e. by introducing network-related
quantities $\rho^{(g)}$ that serve as individual density coefficients.
Specifically, this means modifying formulation (2) by employing
$\rho^{(g)}w^{(g)}(\cdot,\cdot)$, where
$\hat{\rho}^{(g)}=[N^{(g)}(N^{(g)}-1)]^{-1}\sum_{i,j}y_{ij}^{(g)}$ could serve
as an estimate that is independent of the rest of the structure. Such a
framework consequentially might lead to a more balanced comparison strategy.

Beyond the applications presented in the previous section, which all refer to
the situation with two networks, the method could easily be extended to cases
with multiple or even a single network. For example, in the one-sample
setting, to test whether a given network follows a hypothetical distribution
$\mathbb{P}(\boldsymbol{Y}=\cdot\,;\,w(\cdot,\cdot),N)$, we could first align
the network with the theoretical graphon. That is, applying the E-step based
on $w(\cdot,\cdot)$. Given this alignment, we could then turn distributional
assumption (9) into a binomial distribution with the rectangle-specific mean
over $w(\cdot,\cdot)$ as success probability, and, based on that, calculate
the test statistic as in (10).

Besides the testing aspect, our approach could further be used to uncover
relevant differences between networks on the microscopic scale. To be precise,
determining the nodes or edges (present or absent) that contribute most to a
quantified structural discrepancy between networks is possibly interesting in
many situations. This is further elaborated in Section LABEL:sec:diffMicr of
the Supplementary Material.

As a last application case, the joint graphon estimation could further be used
to predict edges between separated networks by considering the cross-sample
probabilities $w^{\text{joint}}(u_{i_{1}}^{(1)},u_{i_{2}}^{(2)})$. This might
be of interest when two (or more) networks are assumed to be samples of a
larger global network. To the best of our knowledge, this has not been pursued
by any other approach so far and hence constitutes a novel perspective. As a
particular hurdle in this framework, the sampling strategy that is supposed
for the drawing of subnetworks needs to be taken into account in the
estimation routine. As far as this is not the Simple Induced Subgraph
Sampling, where one selects a simple random sample of nodes within which all
edges are observed, further adaptations are required. Hence, this lies beyond
the scope of the paper.

SUPPLEMENTARY MATERIAL

Supplementary Manuscript:

(i) Description for deriving differences between networks at the microscopic
level. As exemplary application, the two brain networks from Section 6.2.2 are
considered.
(ii) Replication of the test on functional coactivation networks for two
subgroups of the typical-development group (confer Section 6.2.2).

Python-package for testing on structural equivalence:

Python-package containing the code to perform the comparison methods described
in the paper. The package also contains the preprocessed data of the human
brain functional coactivation networks (see Section 6.2.2). (GNU zipped tar
file)

## References

  * ABIDE I (2013) ABIDE I (2013).  Autism Brain Imaging Data Exchange I. Available online: http://fcon_1000.projects.nitrc.org/indi/abide/ [accessed 06-21-2022].
  * Ali et al. (2014) Ali, W., T. Rito, G. Reinert, F. Sun, and C. M. Deane (2014).  Alignment-free protein interaction network comparison.  Bioinformatics 30(17), i430–i437.
  * Andersen et al. (2016) Andersen, M., J. Dahl, and L. Vandenberghe (2016).  CvxOpt: Open source software for convex optimization (Python). Available online: https://cvxopt.org [accessed 11-03-2022]. Version 1.2.7.
  * Bassett et al. (2018) Bassett, D. S., P. Zurn, and J. I. Gold (2018).  On the nature and use of models in network neuroscience.  Nature Reviews Neuroscience 19(9), 566–578.
  * Bhattacharya et al. (2008) Bhattacharya, K., G. Mukherjee, J. Saramäki, K. Kaski, and S. S. Manna (2008).  The international trade network: weighted network analysis and modelling.  Journal of Statistical Mechanics: Theory and Experiment 2008(02), P02002.
  * Bickel and Chen (2009) Bickel, P. J. and A. Chen (2009).  A nonparametric view of network models and Newman–Girvan and other modularities.  Proceedings of the National Academy of Sciences 106(50), 21068–21073.
  * Bickel et al. (2011) Bickel, P. J., A. Chen, and E. Levina (2011).  The method of moments and degree distributions for network models.  The Annals of Statistics 39(5), 2280–2301.
  * Borgs et al. (2008) Borgs, C., J. T. Chayes, L. Lovász, V. T. Sós, and K. Vesztergombi (2008).  Convergent sequences of dense graphs I: Subgraph frequencies, metric properties and testing.  Advances in Mathematics 219(6), 1801–1851.
  * Burnham and Anderson (2002) Burnham, K. and D. Anderson (2002).  Model selection and multimodel inference: A practical information-theoretic approach. (2nd ed.).  New York: Springer.
  * Butts (2008) Butts, C. T. (2008).  Social network analysis: A methodological introduction.  Asian Journal of Social Psychology 11(1), 13–41.
  * Chan and Airoldi (2014) Chan, S. H. and E. M. Airoldi (2014).  A consistent histogram estimator for exchangeable graph models.  In 31st International Conference on Machine Learning, ICML 2014, pp. 208–216.
  * Craddock et al. (2013) Craddock, C., Y. Benhajali, C. Chu, F. Chouinard, A. Evans, A. Jakab, B. S. Khundrakpam, J. D. Lewis, Q. Li, M. Milham, et al. (2013).  The neuro bureau preprocessing initiative: open sharing of preprocessed neuroimaging data and derivatives.  Frontiers in Neuroinformatics 7, 27.
  * Crossley et al. (2013) Crossley, N. A., A. Mechelli, P. E. Vértes, T. T. Winton-Brown, A. X. Patel, C. E. Ginestet, P. McGuire, and E. T. Bullmore (2013).  Cognitive relevance of the community structure of the human brain functional coactivation network.  Proceedings of the National Academy of Sciences 110(28), 11583–11588.
  * Di Martino et al. (2014) Di Martino, A., C.-G. Yan, Q. Li, E. Denio, F. X. Castellanos, K. Alaerts, J. S. Anderson, M. Assaf, S. Y. Bookheimer, M. Dapretto, et al. (2014).  The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism.  Molecular Psychiatry 19(6), 659–667.
  * Diaconis and Janson (2008) Diaconis, P. and S. Janson (2008).  Graph limits and exchangeable random graphs.  Rendiconti di Matematica e delle sue Applicazioni 28, 33–61.
  * Eagle et al. (2009) Eagle, N., A. Pentland, and D. Lazer (2009).  Inferring friendship network structure by using mobile phone data.  Proceedings of the national academy of sciences 106(36), 15274–15278.
  * Eilers and Marx (1996) Eilers, P. H. and B. D. Marx (1996).  Flexible smoothing with B-splines and penalties.  Statistical Science 11(2), 89–102.
  * Emmert-Streib et al. (2016) Emmert-Streib, F., M. Dehmer, and Y. Shi (2016).  Fifty years of graph matching, network alignment and network comparison.  Information Sciences 346-347, 180–197.
  * Faisal et al. (2017) Faisal, F. E., K. Newaz, J. L. Chaney, J. Li, S. J. Emrich, P. L. Clark, and T. Milenković (2017).  Grafene: Graphlet-based alignment-free network approach integrates 3d structural and sequence (residue order) data to improve protein structural comparison.  Scientific Reports 7(1), 1–15.
  * Gera et al. (2018) Gera, R., L. Alonso, B. Crawford, J. House, J. A. Mendez-Bermudez, T. Knuth, and R. Miller (2018).  Identifying network structure similarity using spectral graph theory.  Applied Network Science 3(1), 1–15.
  * Hoff et al. (2002) Hoff, P. D., A. E. Raftery, and M. S. Handcock (2002).  Latent space approaches to social network analysis.  Journal of the American Statistical Association 97(460), 1090–1098.
  * Holland et al. (1983) Holland, P. W., K. Laskey, and S. Leinhardt (1983).  Stochastic blockmodels: First steps.  Social Networks 5(2), 109–137.
  * Hurvich and Tsai (1989) Hurvich, C. M. and C. L. Tsai (1989).  Regression and time series model selection in small samples.  Biometrika 76(2), 297–307.
  * Kauermann et al. (2013) Kauermann, G., C. Schellhase, and D. Ruppert (2013).  Flexible copula density estimation with penalized hierarchical B-splines.  Scandinavian Journal of Statistics 40(4), 685–705.
  * Kuchaiev et al. (2010) Kuchaiev, O., T. Milenković, V. Memišević, W. Hayes, and N. Pržulj (2010).  Topological network alignment uncovers biological function and phylogeny.  Journal of the Royal Society Interface 7(50), 1341–1354.
  * Latouche and Robin (2016) Latouche, P. and S. Robin (2016).  Variational Bayes model averaging for graphon functions and motif frequencies inference in W-graph models.  Statistics and Computing 26(6), 1173–1185.
  * Leskovec and Krevl (2014) Leskovec, J. and A. Krevl (2014).  SNAP Datasets: Stanford large network dataset collection. Available online: http://snap.stanford.edu/data [accessed 05-27-2022].
  * Leskovec and McAuley (2012) Leskovec, J. and J. McAuley (2012).  Learning to discover social circles in ego networks.  In Advances in Neural Information Processing Systems (NIPS), Volume 25, pp. 539–547.
  * Lovász and Szegedy (2006) Lovász, L. and B. Szegedy (2006).  Limits of dense graph sequences.  Journal of Combinatorial Theory. Series B 96(6), 933–957.
  * Marron and Alonso (2014) Marron, J. S. and A. M. Alonso (2014).  Overview of object oriented data analysis.  Biometrical Journal 56(5), 732–753.
  * Navarro and Segarra (2022) Navarro, M. and S. Segarra (2022).  Joint network topology inference via a shared graphon model.  arXiv preprint arXiv:2209.08223.
  * Newman (2018) Newman, M. E. (2018).  Networks (2nd ed.).  Oxford: Oxford University Press.
  * Nowicki and Snijders (2001) Nowicki, K. and T. A. Snijders (2001).  Estimation and prediction for stochastic blockstructures.  Journal of the American Statistical Association 96(455), 1077–1087.
  * Olhede and Wolfe (2014) Olhede, S. C. and P. J. Wolfe (2014).  Network histograms and universality of blockmodel approximation.  Proceedings of the National Academy of Sciences 111(41), 14722–14727.
  * Onnela et al. (2012) Onnela, J.-P., D. J. Fenn, S. Reid, M. A. Porter, P. J. Mucha, M. D. Fricker, and N. S. Jones (2012).  Taxonomies of networks from community structure.  Physical Review E 86(3), 036104.
  * Pascual-Belda et al. (2018) Pascual-Belda, A., A. Díaz-Parra, and D. Moratal (2018).  Evaluating functional connectivity alterations in autism spectrum disorder using network-based statistics.  Diagnostics 8(3), 51.
  * PCP (2015) PCP (2015).  Preprocessed Connectomes Project: ABIDE. Available online: http://preprocessed-connectomes-project.org/abide/ [accessed 06-21-2022].
  * Pržulj (2007) Pržulj, N. (2007).  Biological network comparison using graphlet degree distribution.  Bioinformatics 23(2), e177–e183.
  * Pržulj et al. (2004) Pržulj, N., D. G. Corneil, and I. Jurisica (2004).  Modeling interactome: Scale-free or geometric?  Bioinformatics 20(18), 3508–3515.
  * Robins et al. (2007) Robins, G., P. Pattison, Y. Kalish, and D. Lusher (2007).  An introduction to exponential random graph (p*) models for social networks.  Social networks 29(2), 173–191.
  * Ruppert et al. (2003) Ruppert, D., M. P. Wand, and R. J. Carroll (2003).  Semiparametric Regression.  Cambridge: Cambridge University Press.
  * Ruppert et al. (2009) Ruppert, D., M. P. Wand, and R. J. Carroll (2009).  Semiparametric regression during 2003–-2007.  Electronic Journal of Statistics 3, 1193–1256.
  * Schwikowski et al. (2000) Schwikowski, B., P. Uetz, and S. Fields (2000).  A network of protein–protein interactions in yeast.  Nature biotechnology 18(12), 1257–1261.
  * Sischka and Kauermann (2022a) Sischka, B. and G. Kauermann (2022a).  EM-based smooth graphon estimation using MCMC and spline-based approaches.  Social Networks 68, 279–295.
  * Sischka and Kauermann (2022b) Sischka, B. and G. Kauermann (2022b).  Stochastic block smooth graphon model.  arXiv preprint arXiv:2203.13304.
  * Snijders and Nowicki (1997) Snijders, T. A. and K. Nowicki (1997).  Estimation and prediction for stochastic blockmodels for graphs with latent block structure.  Journal of Classification 14(1), 75–100.
  * Song et al. (2019) Song, Y., T. M. Epalle, and H. Lu (2019).  Characterizing and predicting autism spectrum disorder by performing resting-state functional network community pattern analysis.  Frontiers in Human Neuroscience 13, 203.
  * Soundarajan et al. (2014) Soundarajan, S., T. Eliassi-Rad, and B. Gallagher (2014).  A guide to selecting a network similarity method.  In SIAM International Conference on Data Mining 2014, SDM 2014, Volume 2, pp. 1037–1045.
  * Subbaraju et al. (2017) Subbaraju, V., M. B. Suresh, S. Sundaram, and S. Narasimhan (2017).  Identifying differences in brain activities and an accurate detection of autism spectrum disorder using resting state functional-magnetic resonance imaging: A spatial filtering approach.  Medical Image Analysis 35, 375–389.
  * Tantardini et al. (2019) Tantardini, M., F. Ieva, L. Tajoli, and C. Piccardi (2019).  Comparing methods for comparing networks.  Scientific Reports 9(1), 1–19.
  * Turlach and Weingessel (2013) Turlach, B. A. and A. Weingessel (2013).  quadprog: Functions to solve quadratic programming problems (R). Available online: https://CRAN.R-project.org/package=quadprog [accessed 11-03-2022]. Version 1.5-5.
  * Tzourio-Mazoyer et al. (2002) Tzourio-Mazoyer, N., B. Landeau, D. Papathanassiou, F. Crivello, O. Etard, N. Delcroix, B. Mazoyer, and M. Joliot (2002).  Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain.  Neuroimage 15(1), 273–289.
  * Wilson and Zhu (2008) Wilson, R. C. and P. Zhu (2008).  A study of graph spectra for comparing graphs and trees.  Pattern Recognition 41(9), 2833–2841.
  * Wolfe and Olhede (2013) Wolfe, P. J. and S. C. Olhede (2013).  Nonparametric graphon estimation.  arXiv preprint arXiv:1309.5936.
  * Wood (2017) Wood, S. N. (2017).  Generalized Additive Models: An Introduction with R.  Boca Raton: CRC Press.
  * Xu et al. (2015) Xu, T., Z. Yang, L. Jiang, X.-X. Xing, and X.-N. Zuo (2015).  A connectome computation system for discovery science of brain.  Science Bulletin 60(1), 86–95.
  * Yaveroğlu et al. (2014) Yaveroğlu, Ö. N., N. Malod-Dognin, D. Davis, Z. Levnajic, V. Janjic, R. Karapandza, A. Stojmirovic, and N. Pržulj (2014).  Revealing the hidden language of complex networks.  Scientific Reports 4(1), 1–9.
  * Yaveroğlu et al. (2015) Yaveroğlu, Ö. N., T. Milenković, and N. Pržulj (2015).  Proper evaluation of alignment-free network comparison methods.  Bioinformatics 31(16), 2697–2704.

## Appendix

###  A Implementation of the Gibbs Sampler

In an iterative joint graphon estimation procedure, the joint posterior
distribution of the node positions given $w^{\text{joint}}(\cdot,\cdot)$ can
be simulated by constructing a Gibbs sampler. We stress that the node
positions are independent between networks and thus the Gibbs sampling
procedure can be conducted for each network separately. The MCMC framework is
then build upon full conditional distribution (5) and can be formulated as
follows. For the successive updating procedure, we consider
$\boldsymbol{u}^{(g),\,<t>}=(u_{1}^{(g),\,<t>},\ldots,u_{N^{(g)}}^{(g),\,<t>})$
to be the current state of the Markov chain. In the $(t+1)$-th step, component
$i$ is then updated according to (5), where all other components remain
unchanged, i.e. $u_{j}^{(g),\,<t+1>}:=u_{j}^{(g),\,<t>}$ for $j\neq i$. To do
so, we propose a new position $u_{i}^{(g),\,*}$ by drawing from a normal
distribution under the application of a logit link. To be precise, we first
calculate

$\displaystyle
v_{i}^{(g),\,<t>}=\operatorname{logit}(u_{i}^{(g),\,<t>})=\log\left(\frac{u_{i}^{(g),\,<t>}}{1-u_{i}^{(g),\,<t>}}\right),$

then we add a normal term in the form of
$v_{i}^{(g),\,*}=v_{i}^{(g),\,<t>}+\operatorname{Normal}(0,\sigma_{v}^{2})$,
and finally we accomplish the retransformation through

$\displaystyle
u_{i}^{(g),\,*}=\operatorname{logit}^{-1}(v_{i}^{(g),\,*})=\frac{\exp(v_{i}^{(g),\,*})}{1+\exp(v_{i}^{(g),\,*})}.$

In this setting, the variance $\sigma_{v}^{2}$ should be chosen such that a
balance between a wide-ranging exploration and a high acceptance rate is
achieved. Given these formulations, the proposal density can be written as

$\displaystyle q(u_{i}^{(g),\,*}|u^{(g),\,<t>}_{i})=$
$\displaystyle\;\frac{\partial u^{(g),\,*}_{i}}{\partial
z^{(g),\,*}_{i}}\phi(z^{(g),\,*}_{i}|z^{(g),\,<t>}_{i})$
$\displaystyle\propto$
$\displaystyle\;\frac{1}{u^{(g),\,*}_{i}(1-u_{i}^{(g),\,*})}$
$\displaystyle\;\cdot\exp\left(-\frac{1}{2}\frac{(\text{logit
}(u_{i}^{(g),\,*})-\text{ logit
}(u^{(g),\,<t>}_{i}))^{2}}{\sigma^{2}}\right),$

which leads to a proposal ratio of

$\displaystyle\frac{q(u^{(g),\,<t>}_{i}|u^{(g),\,*}_{i})}{q(u_{i}^{(g),\,*}|u^{(g),\,<t>}_{i})}=\frac{u^{(g),\,*}_{i}(1-u^{(g),\,*}_{i})}{u^{(g),\,<t>}_{i}(1-u^{(g),\,<t>}_{i})}.$

In combination with the likelihood ration, the acceptance probability of the
proposal, i.e. the probability for setting
$u_{i}^{(g),\,<t+1>}:=u_{i}^{(g),\,*}$, can be calculated through

$\displaystyle\min\left\\{1,\quad\prod_{j\neq
i}\left[\vphantom{\left(\frac{1-w(u^{(g),\,*}_{i},u^{(g),\,<t>}_{j})}{1-w(u^{(g),\,<t>}_{i},u^{(g),\,<t>}_{j})}\right)^{1-y_{ij}}}\right.\right.$
$\displaystyle\left.\left.\left(\frac{w(u_{i}^{(g),\,*},u^{(g),\,<t>}_{j})}{w(u^{(g),\,<t>}_{i},u^{(g),\,<t>}_{j})}\right)^{y_{ij}}\right.\right.$
$\displaystyle\cdot$
$\displaystyle\left.\left.\left(\frac{1-w(u^{(g),\,*}_{i},u^{(g),\,<t>}_{j})}{1-w(u^{(g),\,<t>}_{i},u^{(g),\,<t>}_{j})}\right)^{1-y_{ij}}\right]\frac{u^{(g),\,*}_{i}(1-u^{(g),\,*}_{i})}{u^{(g),\,<t>}_{i}(1-u^{(g),\,<t>}_{i})}\right\\}.$

In case the decision yields a rejection of the proposal, we set
$u_{i}^{(g),\,<t+1>}:=u_{i}^{(g),\,<t>}$. Applying this updating strategy,
which comprises the proposal of a new position plus the decision about its
acceptance, to all $i=1,\ldots,N^{(g)}$ completes one global update. Finally,
we achieve a proper Gibbs sampling routine through consecutively repeating
this global updating scheme. After cutting the burn-in period and applying an
appropriate thinning, this approach yields a sample of the desired joint
posterior distribution of the node positions.

###  B Derivative and Penalization of the B-Spline Function

As has been show in Section 4.2, the log-likelihood of a B-spline function can
be straightforwardly extended towards the situation with multiple datasets.
Given the formulation from (6), the score function can be calculated as

$\displaystyle\boldsymbol{s}(\boldsymbol{\theta})$
$\displaystyle=\left[\frac{\partial\ell(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\right]^{\top}$
$\displaystyle=\sum\limits_{g}\sum\limits_{\begin{subarray}{c}i,j\\\ j\neq
i\end{subarray}}[\boldsymbol{B}_{ij}^{(g)}]^{\top}\left(\frac{y_{ij}^{(g)}}{w_{\boldsymbol{\theta}}^{\text{joint}}(\hat{u}_{i}^{(g)},\hat{u}_{j}^{(g)})}-\frac{1-y_{ij}^{(g)}}{1-w_{\boldsymbol{\theta}}^{\text{joint}}(\hat{u}_{i}^{(g)},\hat{u}_{j}^{(g)})}\right).$

This, in turn, leads to the Fisher information in the form of

$\displaystyle\boldsymbol{F}(\boldsymbol{\theta})$
$\displaystyle=-\mathbb{E}\left(\frac{\partial\boldsymbol{s}(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\right)$
$\displaystyle=\sum\limits_{g}\sum\limits_{\begin{subarray}{c}i,j\\\ j\neq
i\end{subarray}}[\boldsymbol{B}_{ij}^{(g)}]^{\top}\boldsymbol{B}_{ij}^{(g)}\left[w_{\boldsymbol{\theta}}^{\text{joint}}\left(\hat{u}_{i}^{(g)},\hat{u}_{j}^{(g)}\right)\cdot\left(1-w_{\boldsymbol{\theta}}^{\text{joint}}\left(\hat{u}_{i}^{(g)},\hat{u}_{j}^{(g)}\right)\right)\right]^{-1}.$

These results can then be used to implement the Fisher scoring procedure,
where in (7) we additionally add a penalization term to guarantee a smooth
estimation result. For penalizing “neighborhood” elements of the parameter
vector
$\boldsymbol{\theta}=\left(\theta_{11},\ldots,\theta_{1L},\theta_{21},\ldots,\theta_{LL}\right)^{\top}$,
the penalization matrix can be formulated through

$\boldsymbol{P}=\left(\boldsymbol{J}_{L}\otimes\boldsymbol{I}_{L}\right)^{\top}\left(\boldsymbol{J}_{L}\otimes\boldsymbol{I}_{L}\right)+\left(\boldsymbol{I}_{L}\otimes\boldsymbol{J}_{L}\right)^{\top}\left(\boldsymbol{I}_{L}\otimes\boldsymbol{J}_{L}\right),$

where $\boldsymbol{I}_{L}$ is the identity matrix of size $L$ and

$\boldsymbol{J}_{L}=\begin{pmatrix}1&-1&\phantom{-}0&\lx@intercol\hfil\phantom{-}\cdots\hfil\lx@intercol&\phantom{-}0\\\
0&\phantom{-}1&-1&\lx@intercol\hfil\phantom{-}\cdots\hfil\lx@intercol&\phantom{-}0\\\
\vdots&\lx@intercol\hfil\ddots\hfil\lx@intercol&&\phantom{-}\vdots\\\
0&\lx@intercol\hfil\phantom{-}\cdots\hfil\lx@intercol&\phantom{-}0&\phantom{-}1&-1\\\
\end{pmatrix}\in\mathbb{R}^{(L-1)\times L}.$

###  C Choosing the Number and Extent of Rectangles

In order to appropriately test null hypothesis (4), in Section 5, we have
developed an approach that relies on the partition of the graphon’s domain.
According to formulation (8), that involves the number of rectangles, $K$, as
well as their concrete specification in the form of
$[a_{k-1},a_{k})\times[a_{l-1},a_{l})$. In this regard, we emphasize that two
aspects need to be observed. On the one hand, the joint graphon should be
approximately constant within rectangles, requiring
$[a_{k-1},a_{k})\times[a_{l-1},a_{l})$ to be not too expansive. On the other
hand, the amount of edge variables per network falling into these blocks
should be high, which needs rather broad rectangles. Thus, a trade-off between
these two opposed requirements should be reached. In general, we choose $K$ to
grow more slowly than both network dimensions, e.g. scaling as
$\min_{g}\sqrt{N^{(g)}}$. Note that choosing $K=1$ would imply to test whether
the two networks possess the same global density under the assumption of a
joint Erdős–Rényi model. Having determined a suitable value for $K$, we then
simply specify the boundaries of the rectangles through $a_{k}=k/K$ for
$k=0,\ldots,K$. In combination with the subsequent adjustment of the latent
positions as described in Section 4.1, which leads to equidistance of the
estimates $\hat{u}_{i}^{(g)}$, a general lower bound for the amount of
contained nodes per interval, $N_{k}^{(g)}$, can be derived. To be precise, we
can formulate

$N_{k}^{(g)}=\left|\left\\{i\in\\{1,\ldots,N^{(g)}\\}:\frac{i}{N^{(g)}+1}\in[a_{k-1},a_{k})\right\\}\right|\geq\left\lfloor\frac{1}{K}(N^{(g)}+1)\right\rfloor,$

where $\lfloor x\rfloor$ returns the largest integer smaller than or equal to
$x$. Given that, $K$ could also be chosen such that, per network, a prescribed
minimum amount of edge variables per rectangle ($N_{k}^{(g)}N_{l}^{(g)}$ for
$l>k$ and $N_{k}^{(g)}(N_{k}^{(g)}-1)/2$ for $l=k$) is guaranteed.

As a final remark, we emphasize that with regard to the rectangle-based test
statistic, it seems natural to alternatively apply a histogram estimator (Chan
and Airoldi, 2014 or Olhede and Wolfe, 2014). However, the smooth graphon
estimation adapted from Sischka and Kauermann (2022b) considers a global node
ordering which refers not only to separated intervals but to the entire domain
of $[0,1]$. This consequently facilitates the iterative estimation procedure
and thus leads to a more plausible and faster converging node positioning.

###  D Acquiring and Processing of Brain Functional Activation Data

The data we use for analyzing differences in the functional brain activation
are originally provided by the Autism Brain Imaging Data Exchange project
(ABIDE I, 2013, Di Martino et al., 2014). However, we make use of preprocessed
data that are directly accessible through the Preprocessed Connectomes Project
platform (PCP, 2015, Craddock et al., 2013). To be precise, we here apply the
Connectome Computation System pipeline (Xu et al., 2015) and the reduction to
the Automated Anatomical Labeling atlas (Tzourio-Mazoyer et al., 2002). For
each participant, this yields a dataset that consists of activity measurements
over time for $116$ prespecified brain regions (a.k.a. regions of interest).
Based on these temporal activity measurements, we calculate Pearson’s
correlation coefficient between all pairs of brain regions which, per
participant, leads to the corresponding functional connectivity matrix (Song
et al., 2019, Subbaraju et al., 2017). For this analysis, we rely on the data
from New York University, comprising $73$ ASD patients and $98$ TD subjects.
For aggregating these connectivity patterns per clinical group, we apply
Fisher’s transformation to the pairwise correlation coefficients, calculate
their mean for all pairs of brain regions, and finally retransform these means
(Pascual-Belda et al., 2018). This yields for both diagnostic groups a
$116\times 116$ weighted connectivity matrix which we binarize by employing a
threshold of $0.4$. Based on that, the two final networks we obtain both
possess a global density of about $30\%$. With regard to the choice of the
threshold, Song et al. (2019) have found that this has only minor effects when
comparing the networks.