EM-Based Smooth Graphon Estimation Using Bayesian and Spline-Based Approaches

TL;DR

This paper introduces an EM-based method for smooth graphon estimation that combines Bayesian inference and spline regression, allowing for flexible node ordering and uncertainty quantification in network data analysis.

Contribution

It presents a novel EM algorithm integrating Bayesian and spline approaches for nonparametric, smooth graphon estimation without requiring monotonicity constraints.

Findings

Effective in capturing smooth graphons from network data

Allows comparison of different node ordering strategies

Demonstrates robustness through simulations and examples

Abstract

This paper proposes the estimation of a smooth graphon model for network data analysis using principles of the EM algorithm. The approach considers both variability with respect to ordering the nodes of a network and smooth estimation of the graphon by nonparametric regression. To do so, (linear) B-splines are used, which allow for smooth estimation of the graphon, conditional on the node ordering. This provides the M-step. The true ordering of the nodes arising from the graphon model remains unobserved and Bayesian ideas are employed to obtain posterior samples given the network data. This yields the E-step. Combining both steps gives an EM-based approach for smooth graphon estimation. Unlike common other methods, this procedure does not require the restriction of a monotonic marginal function. The proposed graphon estimate allows to explore node-ordering strategies and therefore to compare the common degree-based node ranking with the ordering conditional on the network. Variability and uncertainty are taken into account using MCMC techniques. Examples and simulation studies support the applicability of the approach.

Figures

Tables

Table 1. Table 1: Exemplary graphons considered for simulations. For the second graphon specification Φ ​ ( ) Φ \Phi() and F 𝒩 ​ ( 0 , 0.25 ) subscript 𝐹 𝒩 0 0.25 F_{\mathcal{N}(0,0.25)} denote the cdf’s of the normal distributions with parametrization ( μ , σ 2 ) = ( 0 , 1 ) 𝜇 superscript 𝜎 2 0 1 (\mu,\sigma^{2})=(0,1) and ( 0 , 0.25 ) 0 0.25 (0,0.25) , respectively.

ID	Graphon
1	$\begin{array}{c}{w}_1(u,v)=0.8\left(1-u\right)\left(1-v\right)+0.85\left(u\cdot v\right)\end{array}$
2

Equations

$$U_1,\dots ,U_N\stackrel{\text{i.i.d.}}{\sim }\text{mbox}Uniform[0,1]$$

$$Y_{ij}|U_i=u_i,U_j=u_j\sim \text{mbox}Binomial(1,w(u_i,u_j)).$$

$$g(u)=\int w(u,v)\mathrm{d}v$$

$$U_i\stackrel{\text{i.i.d.}}{\sim }\text{ Uniform}[0,1],i=1,\dots ,N$$

$$\text{mbox}P(Y_{ij}=1|U_i=u_i,U_j=u_j)=w(u_i,u_j)$$

$$f(\mathbf{u}|\mathbf{y})\propto \prod _{\begin{array}{c}i,j\\ j>i\end{array}}w(u_i,u_j{)}^{y_{ij}}(1-w(u_i,u_j){)}^{1-y_{ij}}.$$

$$f_k(u_k|\mathbf{y})\propto \int \dots \int \underset{j>i}{\prod _{i,j}}w(u_i,u_j{)}^{y_{ij}}(1-w(u_i,u_j){)}^{1-y_{ij}}\mathrm{d}{u}_1\dots \mathrm{d}{u}_{k-1}\mathrm{d}{u}_{k+1}\dots \mathrm{d}{u}_N.$$

$$U_{\psi (i)}\le U_{\psi (j)}$$

$$\text{mbox}degree(\hat{\psi }(i))\le \text{mbox}degree(\hat{\psi }(j))$$

$${\hat{u}}_j^{(0)}=\frac{\text{mbox}rank(\text{mbox}degree(j))}{N+1},$$

$${\hat{w}}^{(0)}(u,v)=y_{\hat{\psi }(\lceil uN\rceil )\hat{\psi }(\lceil vN\rceil )},$$

$$f_k(u_k|u_1,\dots ,u_{k-1},u_{k+1},\dots ,u_N,\mathbf{y})\propto \prod _{j\ne k}w(u_k,u_j{)}^{y_{kj}}(1-w(u_k,u_j){)}^{1-y_{kj}}.$$

$$q(u_k^{\ast }|u_k^{<t>})$$

$$\propto \frac{1}{u_k^{\ast }(1-u_k^{\ast })}\exp \left(-\frac{1}{2}\frac{(\text{logit }(u_k^{\ast })-\text{ logit }(u_k^{<t>}){)}^2}{{\sigma }^2}\right),$$

$$\frac{q_k(u_k^{<t>}|u_k^{\ast })}{q_k(u_k^{\ast }|u_k^{<t>})}=\frac{u_k^{\ast }(1-u_k^{\ast })}{u_k^{<t>}(1-u_k^{<t>})}.$$

$$\min \{1,\prod _{j\ne k}\left[{\left(\frac{w(u_k^{\ast },u_j^{<t>})}{w(u_k^{<t>},u_j^{<t>})}\right)}^{y_{kj}}{\left(\frac{1-w(u_k^{\ast },u_j^{<t>})}{1-w(u_k^{<t>},u_j^{<t>})}\right)}^{1-y_{kj}}\right]\frac{u_k^{\ast }(1-u_k^{\ast })}{u_k^{<t>}(1-u_k^{<t>})}\}.$$

$${\bar{u}}_k^{(m)}=\frac{1}{n}\sum _{s=1}^n{u}_k^{<s\cdot N\cdot r>},$$

$${\hat{u}}_j^{(m)}=\frac{\text{mbox}rank({\bar{u}}_j^{(m)})}{N+1}.$$

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

$$\mathbf{\theta }={\left({\theta }_{11},\dots ,{\theta }_{1K},{\theta }_{21},\dots ,{\theta }_{K1},\dots ,,{\theta }_{KK}\right)}^{\top }.$$

$$l(\mathbf{\theta })=\sum _i\sum _{j\ne i}\left[y_{ij}\log \left({\mathbf{B}}_{ij}\mathbf{\theta }\right)+\left(1-y_{ij}\right)\log \left(1-{\mathbf{B}}_{ij}\mathbf{\theta }\right)\right],$$

$$\mathbf{s}(\mathbf{\theta })=\sum _i\sum _{j\ne i}{\mathbf{B}}_{ij}^{\top }\left(\frac{y_{ij}}{w_{\mathbf{\theta }}^{\text{spline}}(u_i,u_j)}-\frac{1-y_{ij}}{1-w_{\mathbf{\theta }}^{\text{spline}}(u_i,u_j)}\right).$$

$$\mathbf{F}(\mathbf{\theta })=\sum _i\sum _{j\ne i}{\mathbf{B}}_{ij}^{\top }{\mathbf{B}}_{ij}{\left[w_{\mathbf{\theta }}^{\text{spline}}\left(u_i,u_j\right)\cdot \left(1-w_{\mathbf{\theta }}^{\text{spline}}\left(u_i,u_j\right)\right)\right]}^{-1}.$$

$$g_{\mathbf{\theta }}^{\text{spline}}(u)=\left[\mathbf{B}(u)\otimes \underset{0}{\overset{1}{\int }}\mathbf{B}(v)\mathrm{d}v\right]\mathbf{\theta }.$$

$$\underset{0}{\overset{1}{\int }}\mathbf{B}(v)\mathrm{d}v$$

$$=\underset{=:\mathbf{A}}{\underbrace{\frac{1}{K-1}\left(\frac{1}{2},1,\dots ,1,\frac{1}{2}\right)}}.$$

$$g_{\mathbf{\theta }}^{\text{spline}}\left({\tau }_l\right)-g_{\mathbf{\theta }}^{\text{spline}}\left({\tau }_{l-1}\right)\ge 0\iff \left[\left(\mathbf{B}\left({\tau }_l\right)-\mathbf{B}\left({\tau }_{l-1}\right)\right)\otimes \mathbf{A}\right]\mathbf{\theta }\ge 0$$

$$\mathbf{L}=\left(\begin{array}{ccccc}1& -1& 0& \dots & 0\\ 0& 1& -1& \dots & 0\\ ⋮& \text{lx@intercol}\text{hfil}\ddots \text{hfil}\text{lx@intercol}& & ⋮& \\ 0& \dots & 0& 1& -1\end{array}\right)\in {\mathbb{R}}^{(K-1)\times K}$$

$$l_{\mathbf{P}}(\mathbf{\theta },\lambda )=l(\mathbf{\theta })-\frac{1}{2}\lambda {\mathbf{\theta }}^{\top }\mathbf{P}\mathbf{\theta },$$

$${\mathbf{s}}_{\mathbf{P}}(\mathbf{\theta },\lambda )=\mathbf{s}(\mathbf{\theta })-\lambda \mathbf{P}\mathbf{\theta }$$

$${\mathbf{F}}_{\mathbf{P}}(\mathbf{\theta },\lambda )=\mathbf{F}(\mathbf{\theta })+\lambda \mathbf{P}.$$

Full text

††Benjamin Sischka is Research Assistant, Department of Statistics, Ludwig-Maximilians-Universität München, 80539 München, Germany (E-mail: benjamin.sischka@stat.uni-muenchen.de). Göran Kauermann is Professor, Department of Statistics, Ludwig-Maximilians-Universität München, 80539 München, Germany (E-mail: goeran.kauermann@stat.uni-muenchen.de).

Keywords: Graphon model; EM algorithm; MCMC; Gibbs sampling; B-spline surface; Social network; Political network; Connectome

1 Introduction

The analysis of network data has achieved increasing interest in the last years. Goldenberg et al. (2010), Hunter et al. (2012), Fienberg (2012) and Salter-Townshend et al. (2012), respectively, published survey articles demonstrating the state-of-the-art in the field. We also refer to Kolaczyk (2009), Kolaczyk and Csárdi (2014) and Lusher et al. (2013) for monographs in the field of statistical network data analysis, see also Kolaczyk (2017). The statistical workhorse models for network data are Exponential Random Graph Models (ERGM), Stochastic Block Models (SBM) and Latent Space Models. ERGMs make use of an exponential family distribution to model the network’s adjacency matrix as a random matrix. This model class was proposed by Frank and Strauss (1986) and is extensively discussed in Snijders et al. (2006). SBMs as well as the latent space models are both model classes which make use of latent quantities relating to the nodes, that are employed to explore the network structure. They were proposed in their nowadays form by Holland et al. (1983) and Hoff et al. (2002), respectively.

A different modeling strategy, which is also based on latent quantities, results through comprehending the network adjacency matrix $\boldsymbol{Y}\in\{0,1\}^{N\times N}$ to be generated by a so called graphon. The graphon as data generating model comes into play by assuming that we draw $N$ random variables

[TABLE]

and, given $U_{i}$ and $U_{j}$, simulate the network entries $Y_{ij}$ conditionally independently through

[TABLE]

The function $w(,)$ is thereby called a graphon (= graph function). In case of undirected networks we additionally require symmetry so that $Y_{ij}=Y_{ji}$ (and hence in principle we assume $w(u,v)=w(v,u)$). We remain with the commonly used convention and exclude so called self-loops, which means $Y_{ii}=0$ for $i=1,\ldots,N$. Apparently, the graphon function $w(,)$ is not unique with regard to the data generating process (2). This is because for any measure-preserving bijection $\varphi:[0,1]\rightarrow[0,1]$ the permuted graphon $w(\varphi(u),\varphi(v))$ yields the same model. More generally, as has been stated by Diaconis and Janson (2008), two graphons $w(,)$ and $w^{\prime}(,)$ represent the same generating model if and only if there exist two measure preserving mappings – not necessarily bijections – $\tilde{\varphi}$ and $\tilde{\varphi}^{\prime}:[0,1]\rightarrow[0,1]$ such that $w(\tilde{\varphi}(u),\tilde{\varphi}(v))=w^{\prime}(\tilde{\varphi}^{\prime}(u),\tilde{\varphi}^{\prime}(v))$ for almost all $(u,v)\in[0,1]^{2}$. Some papers therefore add a further attribute to achieve uniqueness, e.g. see Bickel and Chen (2009) or Yang et al. (2014). The common setting to do so is to postulate that

[TABLE]

is strictly increasing in $u$, which leads to the so called canonical representation of the graphon $w(,)$. Note that the distribution of $g(U_{i})$ with $U_{i}$ following (1) can be interpreted as (asymptotic) distribution of the degree proportion. However, the additional condition (3) incorporated to tackle the identifiability issue can be problematic, especially in combination with other postulations. For reconstruction reasons commonly some kind of smoothness is assumed in the sense that $w(,)$ satisfies (at least piecewise) some Lipschitz or Hölder conditions of continuity, e.g. see Olhede and Wolfe (2014), Chan and Airoldi (2014) or Gao et al. (2015a). But then, considering an arbitrary continuous graphon, the transformation into its canonical representation neither usually retains the continuity nor are there necessarily any guarantees that such a canonical representation exists. Taking only graphon functions into account which are continuous in their canonical representation would apparently be a strong restriction of the generality of graphon models. We therefore consider the more general model class which does not include the canonical constraint. Still, we also discuss estimation under the restriction (3).

Graphon estimation for modeling network data has recently found attention in the statistical literature. A general discussion is given in Orbanz and Roy (2014) and You (2020) recently launched an R package for graphon estimation. Graphons can be related to ERGMs, at least for simple statistics like two-star or triangles, as shown in Diaconis and Chatterjee (2013). He and Zheng (2015) make use of this connection and propose to use asymptotic properties of graphons to derive estimates in high dimensional ERGMs. Moreover, graphon models can be seen as a generalized model class which also includes SBMs and latent space models (e.g. see Bickel et al., 2011). Wolfe and Olhede (2013) and Yang et al. (2014) discuss non-parametric estimation of graphons including tests on the validity of prespecified graphon shapes, see also Chan and Airoldi (2014), Airoldi et al. (2013) or Olhede and Wolfe (2014). One of the most propagated strategies in the graphon estimation literature is to approximate the graphon through a SBM, see Choi and Wolfe (2014), Choi (2017) or Klopp et al. (2017). Gao et al. (2015a) discuss optimal graphon estimation for the SBM approximation. SBMs are, however, by definition discontinuous models, meaning that the corresponding graphon function is discontinuous and hence not smooth. In this paper, in contrast, we focus on smooth and hence continuous graphon functions. For a general discussion on graphons we refer to Borgs et al. (2008), Lovász (2012), Diaconis and Janson (2008) or Bickel and Chen (2009).

Estimating graphons can generally be reduced to probability matrix estimation, where the goal is to gain information about the specific pairwise probabilities $\mbox{P}(Y_{ij}=1|U_{i}=u_{i},U_{j}=u_{j})=w(u_{i},u_{j})$. From this perspective, the $U_{i}$ are not well-defined components. This reduction – considering the graphon values only at unspecific positions – allows to circumvent the identifiability issue from above. Chatterjee (2015) provides convergence results for general matrix estimation using singular value decomposition, see also Gao et al. (2015a), Zhang et al. (2015) or Xu (2017). Gao et al. (2016) extends this towards partially observed matrices/networks, see also Gao et al. (2015b) for a Bayesian approach. The cited work, however, yields estimates of specific graphon values at unspecific positions and thus does not lead to a smooth graphon estimate on the domain $[0,1]^{2}$ which is focused in this work.

In this paper we propose to use penalized linear B-splines for graphon estimation. This borrows ideas suggested in Kauermann et al. (2013) for copula estimation, since B-splines easily allow to accommodate side constraints such as symmtery in form of $w(u,v)=w(v,u)$ or, if required, condition (3) for the resulting estimate. This, in contrast, is difficult to accommodate in histogram or kernel based estimation. Penalized estimation with B-splines has thereby a long standing tradition in smooth estimation, starting with Eilers and Marx (1996) and Ruppert et al. (2003, 2009), see also Wood (2017). We extend the idea here to graphon estimation. However, before this smoothing approach can be applied, we also need to fill the lack of information about the $U_{i}$. On the other hand, conditional on the observed network data $\boldsymbol{y}=\{y_{ij}\}_{i,j=1,\ldots,N}$, a presumption about the positions of the $U_{i}$ can only be made in relation to $w(,)$. Since both the implicitly supposed parameters of the B-spline approximation of $w(,)$ and the latent $U_{i}$ are unknown, this is a typical task for an EM type algorithm. Moreover, MCMC techniques can be applied in this context to approximate numerically the complex posterior distribution of the $U_{i}$, which together yields a MCEM algorithm.

The paper is organized as follows. Section 2 displays the main ideas of pursuing an EM based algorithm for smooth graphon estimation. Section 3 describes the procedure in detail. Section 4 and Section 5 showcase results for both simulations and real-world data examples. A discussion concludes the paper.

2 Graphon Representation and EM Motivation

We assume that the graphon $w:[0,1]^{2}\rightarrow[0,1]$ is a smooth function, meaning that $w(,)$ is sufficiently differentiable in both arguments. We call $w(,)$ to have a canonical representation if $g(u)=\int w(u,v)\mathop{}\!\mathrm{d}v$ is strictly increasing, see, among others, Chan and Airoldi (2014). We assume further that $w(,)$ is symmetric and generates a network of size $N$ through the following process. For $N$ independent uniform variables

[TABLE]

we obtain the symmetric network through

[TABLE]

for $1\leq i<j\leq N$, where $Y_{ji}=Y_{ij}$ and $Y_{ii}=0$. The variables $U_{i}$ remain unobservable and as data we only obtain the observed network $\boldsymbol{y}$. This implies in particular that no information on the values of $U_{i}$ is given and hence estimation of $w(,)$ occurs to be difficult. We propose to tackle the estimation problem using an EM algorithm. In that sense, we calculate (or rather approximate by simulating) $\mbox{E}(\boldsymbol{U}|\boldsymbol{y})$ from which we derive an appropriate ordering, giving the E-step, which in turn allows to estimate $w(,)$ using smoothing techniques, providing the M-step. For the E-step we take a Bayesian view by looking at the posterior density of $\boldsymbol{U}=(U_{1},\ldots,U_{N})$ given $\boldsymbol{y}$. Note that since $U_{1},\ldots,U_{N}$ are i.i.d. uniform we have

[TABLE]

If we look at the univariate distribution of a single variable $U_{k}$ given the entire network $\boldsymbol{y}$, this results through

[TABLE]

Apparently, both the joint and the marginal posterior distribution are too complex to be calculated analytically, in particular if $N$ is large. We will therefore explore (5) by pursuing a Bayesian approach using MCMC techniques. Details are given in the next sections. Letting $\hat{u}_{k}^{(m)}$ denote the rescaled ranks of the posterior means (see equation (9) below) of the MCMC samples in the $m$th step of the EM algorithm, we then maximize the likelihood resulting from (4) with $u_{i}$ replaced by $\hat{u}_{i}^{(m)}$. To do so, we make use of penalized B-spline smoothing, meaning that we replace the unknown graphon $w(,)$ by an approximate spline base representation $\boldsymbol{B}(,)\boldsymbol{\theta}$, where $\boldsymbol{B}(,)$ is a bivariate spline and $\boldsymbol{\theta}$ is the vector of unknown spline coefficients to be estimated in the M-step. Both steps the E- and the M-step will be introduced in detail below. Before doing so, however, we propose a simple first E-step, which also serves as initial step of the corresponding EM iteration.

Even though the probability model (4) used for the construction of networks is simple, it can not directly be used for estimation. The reason is that variables $U_{i}$ are unobservable and hence can not directly be employed to estimate the graphon $w(,)$. Instead, in the recent literature the graphon $w(,)$ is usually estimated by smoothing the observed adjacency matrix $\boldsymbol{y}$. The preceding rearrangement of the nodes thereby should effect an appropriate reordering of the rows and columns which in turn should reflect the ordering of the $U_{i}$. Let $\psi:\{1,\ldots,N\}\rightarrow\{1,\ldots,N\}$ be a permutation such that

[TABLE]

with $i<j$. This means $U_{\psi(i)}=U_{(i)}$, where $U_{(1)}\leq U_{(2)}\leq\ldots\leq U_{(N)}$ define the ordered variables $U_{i}$. Note that since $U_{i}$, $i=1,\ldots,N$ are not observable, we can also not observe $\psi()$ which therefore needs to be estimated. A common approach is to make use of the degree. Let therefore $\hat{\psi}:\{1,\ldots,N\}\rightarrow\{1,\ldots,N\}$ be a permutation such that

[TABLE]

for $i<j$. Note that $\hat{\psi}()$ can serve as an initial estimate for $\psi()$, although it requires the assumption that $g()$ from (3) is strictly increasing. We therefore use it as initial E-step and define the corresponding resulting initial prediction for $U_{j}$ based on this simple sorting through

[TABLE]

where $\mbox{{rank}}(\mbox{{degree}}(j))$ is the rank from smallest to largest of the $j$th element of the tuple $(\mbox{{degree}}(i):\;\allowbreak i=1,\ldots,N)$. This is equivalent to define $\hat{u}_{\hat{\psi}(j)}^{(0)}=j/(N+1)$, where $i/(N+1)$, $i=1,\ldots,N$ represent the expected values of $N$ ordered independently $\text{Uniform}[0,1]$ distributed variables. Chan and Airoldi (2014) prove – in case of graphons with canonical representation – asymptotic convergence rates for $\hat{\psi}()$ towards $\psi()$, meaning that $|\psi(j)-\hat{\psi}(j)|/N\rightarrow 0$ for all $j=1,\ldots,N$.

We can now replace $w(,)$ in (4) by its empirical version $\hat{w}^{(0)}(,)$ which is defined through

[TABLE]

where $\lceil uN\rceil$ defines the smallest value which is greater or equal to $uN$. Note that $\hat{w}^{(0)}(,)$ just mimics the ordered adjacency matrix scaled towards the unit square. These calculations provide the initial steps in the subsequently introduced EM algorithm.

3 EM Algorithm for smooth Graphons

3.1 Bayesian Approach for the E-Step

We pursue the E-step by exploiting the posterior distribution of $\boldsymbol{U}$. This is done by constructing an appropriate MCMC Gibbs sampling scheme based on the full posterior distribution of $U_{k}$. Note that by conditioning on $\boldsymbol{y}$ and all $U_{j}$ except of $U_{k}$ one gets

[TABLE]

We pretend in this section that the graphon $w(,)$ is known, which is the general setting in the E-step. This allows to easily sample from (8) using Gibbs sampling as MCMC technique. To do so, we assume $\boldsymbol{u}^{<t>}=(u^{<t>}_{1},\ldots,u^{<t>}_{N})$ to be the current state of the Markov chain. To update the $k$th component we then set $u^{<t+1>}_{j}:=u^{<t>}_{j}$ for $j\neq k$, while component $u_{k}$ is updated by drawing from (8). To pursue this, we make use of a normal proposal using a logit link. To be specific, let $z^{<t>}_{k}=\log(u^{<t>}_{k}/(1-u^{<t>}_{k}))=\text{ logit }(u^{<t>}_{k})$. We then propose to draw $z^{*}_{k}=z^{<t>}_{k}+\mbox{Normal}(0,\sigma^{2})$ and set u_{k}^{*}=\text{ logit}^{-1}(z_{k}^{*})=\exp(z^{*}_{k})\big{/}(1+\exp(z_{k}^{*})). Hence, the proposal density for $U_{k}$ is proportional to

[TABLE]

where $\phi()$ is the standard normal density. Consequently, the ratio of proposals equals

[TABLE]

The proposed value $u^{*}_{k}$ is accepted (and hence we set $u^{<t+1>}_{k}:=u_{k}^{*}$) with probability

[TABLE]

If we do not accept $u_{k}^{*}$, we set $u^{<t+1>}_{k}:=u^{<t>}_{k}$. Based on the resulting Markov chain we can estimate the posterior mean $\mbox{E}(U_{k}|\boldsymbol{y})$ by taking the mean of the simulated values, observing an appropriate burn-in phase of the algorithm. To be specific, we use the MCMC sequence to estimate the posterior mean in the $m$th step of the EM algorithm through

[TABLE]

where $r\in\mathbb{N}$ describes a thinning factor and $n$ is the number of MCMC states which are taken into account. We then use the ranks of the posterior means to reorder the network matrix accordingly. This corresponds to setting the missing values of $U_{j}$ according to (7) to

[TABLE]

The Gibbs sampling approach is straightforward and simple but requires the knowledge of the graphon $w(,)$. Apparently, in the EM algorithm we replace $w(,)$ in the formula above by the current estimate resulting through the M-step, which is described subsequently. We define the final estimate resulting through (9) as $\hat{u}_{j}^{\text{{EM}}}$.

3.2 Spline based Graphon Estimation for the M-Step

3.2.1 Linear B-Splines

For smooth estimation of the graphon $w(,)$ we first formulate a spline based approximation through

[TABLE]

where $\boldsymbol{B}(u)\in\mathbb{R}^{1\times K}$ is a linear B-spline basis on $[0,1]$, normalized to have maximum value 1, see Figure 1. Parameter vector $\boldsymbol{\theta}\in\mathbb{R}^{K^{2}\times 1}$ is indexed through

[TABLE]

Using (10) we obtain the likelihood

[TABLE]

where $\boldsymbol{B}_{ij}=\boldsymbol{B}(u_{i})\otimes\boldsymbol{B}(u_{j})$. Taking the derivative leads to the score function

[TABLE]

Moreover, taking the expected second order derivative leads to the Fisher matrix

[TABLE]

Our intention is to maximize $l(\boldsymbol{\theta})$, which could be simply done by Fisher scoring. The resulting maximizer does, however, not by default lead to a proper estimate, meaning to fulfill symmetry and boundedness. Furthermore, in case we aim to estimate a graphon with canonical representation, we need to incorporate the constraint from (3). We therefore impose additional (linear) side constraints on $\boldsymbol{\theta}$. Considering the identifiability issue we get the marginal function with (10) through

[TABLE]

For normalized B-splines we can easily calculate the integral and for a standardized basis with equidistant knots we obtain

[TABLE]

This allows to rewrite (11) to $g_{\boldsymbol{\theta}}^{\text{{spline}}}(u)=\left[\boldsymbol{B}(u)\otimes\boldsymbol{A}\right]\boldsymbol{\theta}$. Hence, the marginal function $g_{\boldsymbol{\theta}}^{\text{{spline}}}()$ is also expressed as a linear B-spline and a monotonicity constraint is easily accommodated by postulating monotonicity at the knots $\tau_{1},\ldots,\tau_{K}$. That is we need

[TABLE]

for $l=2,\ldots,K$, which is a simple linear constraint on the coefficient vector.

Imposing symmetry on the graphon can also easily be accommodated as linear constraints $\theta_{pq}=\theta_{qp}$ for $p\neq q$. Finally, we need $w_{\boldsymbol{\theta}}^{\text{{spline}}}(,)$ being bounded to $[0,1]$, which is again a simple linear constraint. All in all, we can write the side constraints as $\boldsymbol{C}\boldsymbol{\theta}\geq\boldsymbol{0}$ and $\boldsymbol{D}\boldsymbol{\theta}=\boldsymbol{0}$ for matrices $\boldsymbol{C}$ and $\boldsymbol{D}$ chosen accordingly, where the constraint from (12) for the canonical representation can be added if required. With the above linear constraints and the maximization task of $l(\boldsymbol{\theta})$ we obtain an (iterated) quadratic programming problem, which can be solved using standard software (see e.g. Andersen et al., 2004 or Turlach and Weingessel, 2013).

3.2.2 Penalized Estimation

Following ideas from the penalized spline estimation (see Eilers and Marx, 1996 or Ruppert et al., 2009) we may additionally impose a penalty on the coefficients to achieve smoothness. This is necessary since we intend to choose $K$ large and unpenalized estimation will lead to wiggled estimates. We refer to Eilers and Marx (1996) for a motivation of penalized spline estimation. To do so, we penalize the difference between “neighbouring” elements of $\boldsymbol{\theta}$ to achieve smoothness. Let therefore

[TABLE]

be the first order difference matrix. We then penalize $\left[\boldsymbol{L}\otimes\boldsymbol{I}\right]\boldsymbol{\theta}$ and $\left[\boldsymbol{I}\otimes\boldsymbol{L}\right]\boldsymbol{\theta}$, where $\boldsymbol{I}$ is the identity matrix of appropriate size. This is leading to the penalized likelihood

[TABLE]

where $\boldsymbol{P}=\left(\boldsymbol{L}\otimes\boldsymbol{I}\right)^{\top}\left(\boldsymbol{L}\otimes\boldsymbol{I}\right)+\left(\boldsymbol{I}\otimes\boldsymbol{L}\right)^{\top}\left(\boldsymbol{I}\otimes\boldsymbol{L}\right)$ and $\lambda$ serves as smoothing parameter. The corresponding penalization score function is given through

[TABLE]

and the penalized Fisher matrix in the form

[TABLE]

The estimate apparently depends on the penalty parameter $\lambda$, which is suppressed in the notation. Setting $\lambda\rightarrow 0$ gives an unpenalized fit while setting $\lambda\rightarrow\infty$ leads to a constant graphon, i.e. an Erdős-Rényi model. The smoothing parameter $\lambda$ therefore needs to be chosen data driven. We here follow Kauermann et al. (2013) and make use of the Akaike Information Criterion (AIC) (Hurvich and Tsai, 1989, see also Burnham and Anderson, 2010). To do so, we define the corrected AIC through

[TABLE]

where $\hat{\boldsymbol{\theta}}_{\boldsymbol{P}}$ is the penalized parameter estimate and $df(\lambda)$ is the degree of the model, which we define in the usual way as trace of the product of the inverse penalized Fisher matrix and the unpenalized Fisher matrix. To be specific

[TABLE]

We subsequently denote with $\hat{w}^{(1)}(,)$ and $\hat{w}^{\text{{EM}}}(,)$ the penalized B-spline estimates of $w_{\boldsymbol{\theta}}^{\text{{spline}}}(,)$ in the first and the final EM step, meaning that $\hat{w}^{(1)}(,)$ and $\hat{w}^{\text{{EM}}}(,)$ are based on $\hat{\boldsymbol{u}}^{(0)}$ and $\hat{\boldsymbol{u}}^{\text{{EM}}}$, respectively.

As it is well known, the EM algorithm can be trapped at a local maximum of the likelihood (or, in this case, at a local minimum of the corrected AIC). Thus, depending on the specific data situation, it might be recommendable to repeat the algorithm several times to achieve a better fit.

3.3 Information on Ranking

Note that $U_{i}$ is uniform and it is helpful to order $U_{i}$ such that $U_{(1)}\leq U_{(2)}\leq\ldots\leq U_{(N)}$. Considering the degree based ordering $\hat{\psi}()$ from (6) as representative, this allows to define the (degree related) ranking density $f_{\hat{\psi}(k)}(u_{\hat{\psi}(k)}|\boldsymbol{y})$. For the sake of simplicity, we subsequently collapse $\hat{\psi}(k)$ to $(k)$ so that henceforth $f_{(k)}(u_{(k)}|\boldsymbol{y})$ and $U_{(k)}$ refer to the node with the $k$th lowest degree. For other quantities the notation applies accordingly. This ranking density, however, again is difficult or even impossible to calculate analytically. Note that the full posterior density is given through

[TABLE]

where $y_{(k)(j)}$ is the network entry which refers to the nodes $k$ and $j$ after the permutation according to the degree based ordering. In that regard, the MCMC sequence (after appropriate thinning) provides information about the posterior distribution of $\boldsymbol{U}$ given the network $\boldsymbol{y}$. Hence, for the marginal posterior distribution of $U_{(k)}$ we can follow a Monte Carlo integration approach, see e.g. Gelfand and Smith (1990, sec. 2.2 and 2.3), and calculate

[TABLE]

where $\boldsymbol{u}^{<t>}_{(-k)}=(u^{<t>}_{(1)},\ldots,u^{<t>}_{(k-1)},u^{<t>}_{(k+1)},\ldots,u^{<t>}_{(N)})$ is the $t$th state of the Gibbs sampling sequence without the $k$th component after degree related permutation, $r\in\mathbb{N}$ describes a thinning factor and $n$ is the number of MCMC states which are taken into account. For this purpose, the unspecified normalizing constant in (13) can be approximated through a Riemann sum, since we assume $w(,)$ to fulfill certain continuity properties. Considering $w(,)$ as unknown again, we employ $\hat{w}^{(1)}(,)$ or $\hat{w}^{\text{{EM}}}(,)$ as estimate in (13) and obtain $\hat{f}_{(k)}^{(1)}(u_{(k)}|\boldsymbol{y})$ or $\hat{f}_{(k)}^{\text{{EM}}}(u_{(k)}|\boldsymbol{y})$ as estimate of (14), respectively.

4 Simulation Studies

For evaluating our method we at first consider networks generated from a known ground truth. More precisely, for each of the two graphons from Table 1 we simulate networks with dimension $N=500$ using the data generating process (4).

The first graphon has a canonical representation, while the second is more general and does not fulfill that $g()$ in (3) is strictly increasing. We start with Graphon 1 and demonstrate the benefits of iteratively applying the E- and M-step instead of ordering the nodes just with respect to their degree, i.e. based on $\hat{\psi}()$, and applying the M-step only once. We call the latter a “one-step” estimator. This “one-step” estimate $\hat{w}^{(1)}(,)$ is shown in the top right panel in Figure 2.

The fit is apparently not convincing, which we want to make more evident by looking at the ranking information. First, we plot $\hat{u}^{(0)}_{i}$ as defined in (7) against the true values $u_{i}$. There is no concordance visible, which exhibits the problem with a simple ranking of the nodes based on the degree. We run an MCMC sample according to Section 3.1 to obtain information about the ranking, as discussed in Section 3.3. The MCMC sample allows to estimate the posterior distribution of $U_{(k)}$ as proposed in (14), which for three exemplary values of $k$ is plotted in the bottom right plots of Figure 2. The corresponding initial estimates $\hat{u}_{(k)}^{(0)}$ are included as vertical dashed lines and it can be seen that they are not well represented by the posterior distributions. We can conclude that this network, which has been generated from $w_{1}(,)$, is not suitable to sort by degree for achieving an appropriate estimate of the true node ordering with respect to the data generating model. Hence, also the “one-step” estimation approach results in a very poor fit of the underlying true graphon. We therefore iterate between the E- and the M-step.

The converged EM estimate is shown in the top right plot in Figure 3.

Regarding the final EM based estimates $\hat{u}_{i}^{\text{{EM}}}$, the comparison with the true values $u_{i}$ reveals a convincing concordance. We also consider again the posterior distribution of $U_{(k)}$ for the same indices as above (bottom right plots), which indicates a plausible positioning. As a conclusion, the proposed EM algorithm provides appropriate results even if the initial ordering of the nodes by their degree is not adequate. Note that the canonical restriction from (12) has not been included here although $w_{1}(,)$ has canonical representation, since the shown unrestricted estimate turned out to be similarly good.

Finally, we explore the convergence behavior of the algorithm, which is illustrated in Figure 4 for some values of the graphon at exemplary specific positions. We see that after about 8 to 12 steps a reasonable convergence occurs.

For a further evaluation, we now have a look at Graphon 2, which does not provide a canonical representation and thus is completely intractable for simple degree based estimation methods. This simulation seeks to demonstrate that our proposed method can also handle more general graphons. The marginal function $g()$ as defined in (3) is constant at $0.25$, meaning that the degree is fully uninformative with regard to the node ordering. Nonetheless, we here remain with the degree based initialization as proposed above to show that the EM algorithm is capable of appropriately estimating the underlying structure even under bad initialization. We exclude the side constraint from (12), since we now explicitly want to enable flexible marginal functions. The resulting graphon estimate for a simulated network of dimension $N=500$ is depicted in the top right plot in Figure 5.

The structure of the true graphon (top left) is fully captured, which is accompanied by a convincing concordance between the estimated $\hat{u}^{\text{{EM}}}_{i}$ and the true $u_{i}$. In accordance with that, the estimates $\hat{u}^{\text{{EM}}}_{(k)}$ are well covered by the corresponding ranking densities, which is illustrated for three selected indices at the stacked plots at the bottom right. To further showcase the performance of our method we repeat the estimation procedure several times based on the same simulated network but with differing random initialization. More precisely, instead of making use of the degree based ordering from (7), we set $\hat{\boldsymbol{u}}^{(0)}$ as a random permutation of $(i/(N+1):\;\allowbreak i=1,\ldots,N)$. This also exemplary characterizes the appearance of the EM algorithm being trapped at local minima of the AIC. Figure 6 illustrates the estimation results for six such repetitions.

Although we start for each run with a completely uninformative random initialization, in most of the final estimates the structure of the original graphon can be fully recognized instantly. Yet, the estimates in the top right and the bottom left panel reveal a somehow differing structure. Having a more precise look, they exhibit the appearance of segment swaps and reversals with respect to the original graphon representation. For instance, at the graphon estimate in the bottom left panel the upper and the lower part of the domain $[0,1]$ have been swapped and, in addition, the lower part has been reversed. However, as mentioned before, applying permutations to the graphon function has no effect on the model specification itself. Looking at the resulting AIC values beneath the graphon estimates, we nevertheless see that the smallest values result for smooth graphon estimates, while the occurrence of a “jump” due to only piecewise correct reordering of the nodes lead to increased AIC values. But still, the structure of the true graphon model is always well captured, albeit possibly in a different form of representation. Summarizing over all the estimation results, our algorithm yields very good estimates for Graphon 2, although an uninformative initial node ordering has been applied.

The example illustrates that a random reordering of the nodes for the initial E step can be applied in combination with the AIC in order to obtain an optimal smooth graphon estimate.

5 Real-World Data Examples

We complete the paper by analyzing real-world network data. To do so, we consider networks from three different domains, namely from sociology, political science, and neuroscience. For the EM based estimation we use degree ordering as initialization where the structure seems appropriate. Otherwise we apply random initialization and select the best outcome over several repetitions.

5.1 Facebook Network

A very common application and one of the roots of network sciences are friendship networks. We here consider a Facebook ego network which has been collected by McAuley and Leskovec (2012) and is available on the Stanford Large Network Dataset Collection (Leskovec and Krevl, 2014). This ego network with $333$ actors is plotted in Figure 7 with two different colorings.

In the left panel the illustrated ordering is based on the degree, while on the right it refers to the final EM result. Apparently, the ordering given by the final EM result seems much more appropriate with respect to the network structure. Thus, the iterative approach improves the initial degree ordering markedly. Moreover, the inherent structure of the network on the right in Figure 7 can to some extent be recognized instantly in the corresponding final graphon estimate $\hat{w}^{\text{{EM}}}(,)$, which is shown on the left in Figure 8.

For example, the bundle of nodes in the center of the lower network part with roughly $\hat{u}_{i}^{\text{{EM}}}\in[0.65,1]$ and the connectivity pattern among themselves is captured in the graphon estimate in the intense section on the bottom right. In accordance with this, also the estimates $\hat{u}_{(k)}^{\text{{EM}}}$ for some selected indices are adequately represented by the corresponding posterior distributions of $U_{(k)}$, which underlines the appropriateness of the graphon estimate $\hat{w}^{\text{{EM}}}(,)$.

5.2 Military Alliance Network

As second real-world network example, we consider strong military alliances among the world’s nations. For that purpose, we use data from the Alliance Treaty Obligations and Provisions project (Leeds et al., 2002), which provide information about all kind of military alliance agreements over an extensive period. For a suitably reduced network we here define the presence of a strong military alliance when states have entered into an offensive or defensive pact, meaning when they have signed a treaty which forces the one country to intervene by military active support if the other country comes into a conflict with offensive or defensive actions, respectively. Furthermore, we truncate the data to agreements that were in force in 2016 as the most recent available year. The best final estimation results of the EM approach over several repetitions with random initial node ordering are illustrated in Figure 9.

The graphon estimate in the top right panel exhibits a very pronounced assortative structure, meaning that predominantly pairs of nodes are linked whose latent quantities are close. The final node ordering, which is visualized by the coloring in the top left network, appears reasonable with respect to the network structure. Transferring this coloring to the world map reveals a strong conformity between the geographic closeness and the closeness of the latent quantities. Together with the assortative structure, this means that countries which are geographically close form similar military alliances and further are more likely to be allied with each other.

5.3 Human Brain functional Coactivation Network

To conclude the real-world examples, we consider a human brain functional network which has been constructed from a meta-analysis by Crossley et al. (2013) and which is available in the Brain Connectivity Toolbox (Rubinov and Sporns, 2010a, see also Rubinov and Sporns, 2010b for detailed description). Their weighted network matrix represents the “estimated […] similarity (Jaccard index) of the activation patterns across experimental tasks between each pair of $638$ brain regions” (Crossley et al., 2013). To obtain a binary adjacency matrix out of it, we simply apply a threshold of above zero, meaning that we include a link between each pair of brain regions which is at least during one task activated simultaneously. The resulting network with a density of approximately $9.2\%$ is depicted in the top left panel in Figure 10, where the coloring refers to $\hat{u}_{i}^{\text{{EM}}}$.

For the evaluation, we again consider the best outcome of the EM algorithm over several repetitions with random initialization. The graphon estimate in the top right plot also reveals an assortative structure but, in addition, exhibits a conspicuous pattern of functional coactivations of some segments, which are separated with respect to the latent dimension. Regarding the spatial positions of the brain regions, the lower three plots show a strong relation between closeness in anatomical space and closeness of the $\hat{u}_{i}^{\text{{EM}}}$. Nevertheless, there also seem to be areas which have a similar color pattern (indicating a similar coactivation pattern) but are spatially separated. Indeed, the interaction and coactivation of distant brain areas is a familiar phenomena in neuroscience and can also be seen in Crossley et al. (2013), who pursue a different modeling strategy. Overall, we can demonstrate that the EM approach for estimating smooth graphons provides additional insight into network structures.

6 Discussion

The paper proposes a novel estimation routine for graphon estimation which explicitly takes the variability of ordering the nodes into account. The proposed semi-parametric estimation based on (linear) B-splines allows to incorporate relevant properties in the estimation and also uniqueness restrictions if required. The Bayesian approach relying on Gibbs sampling illuminates the uncertainty about the degree and its distribution. Both steps combined give an EM type algorithm which allows for flexible graphon estimation even in large networks. The approach outperforms available routines in two aspects. First, the B-spline estimate can guarantee a smooth and continuous outcome and also a unique representation of the graphon if required. Secondly, based on the Bayesian formulation and the EM algorithm one can assess the amount of uncertainty for ordering the nodes based on their degree or based on any other strategy.

The proposed approach can also be used in other related models like stochastic block models (SBM), where we assume that nodes cluster and form within and between the clusters simple Erdős-Rényi models. Apparently, SBMs do not have a smooth underlying graphon structure so that one needs to adjust and further develop the results in this paper.

Overall, graphon estimation also provides an interesting tool for network visualizations, as demonstrated in the examples. In this respect, it is more than a modeling exercise but may also serve as tool for explanatory network data analysis.

Acknowledgments

The project was partially supported by the European Cooperation in Science and Technology [COST Action CA15109 (COSTNET)]. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. Declarations of interest: none.