# Bayesian classification, anomaly detection, and survival analysis using
network inputs with application to the microbiome

Nathaniel Josephs  Lizhen Lin  Steven Rosenberg  Eric D. Kolaczyk  Department
of Mathematics and Statistics, Boston University  Department of Applied and
Computational Mathematics and Statistics, The University of Notre Dame

###### Abstract

While the study of a single network is well-established, technological
advances now allow for the collection of multiple networks with relative ease.
Increasingly, anywhere from several to thousands of networks can be created
from brain imaging, gene co-expression data, or microbiome measurements. And
these networks, in turn, are being looked to as potentially powerful features
to be used in modeling. However, with networks being non-Euclidean in nature,
how best to incorporate them into standard modeling tasks is not obvious. In
this paper, we propose a Bayesian modeling framework that provides a unified
approach to binary classification, anomaly detection, and survival analysis
with network inputs. We encode the networks in the kernel of a Gaussian
process prior via their pairwise differences and we discuss several choices of
provably positive definite kernel that can be plugged into our models.
Although our methods are widely applicable, we are motivated here in
particular by microbiome research (where network analysis is emerging as the
standard approach for capturing the interconnectedness of microbial taxa
across both time and space) and its potential for reducing preterm delivery
and improving personalization of prenatal care.

Multiple Networks,

Gaussian Process,

Graph Kernel,

Frobenius Distance,

Spectral Distance,

Random Walk Kernel,

Microbiome Networks,

Preterm Delivery,

###### keywords:

, ,
and

##  1 Introduction

Preterm delivery is a global health problem and hence an important area of
research. According to the World Health
Organization,111https://www.who.int/en/news-room/fact-sheets/detail/preterm-
birth it is the leading cause of death among infants worldwide, but it is
estimated that over 75% of preterm infants would survive with appropriate
intervention with a survival rate higher among infants born the latest.
Therefore, it is important to be able to identify women at a high risk of
preterm delivery and also be able to predict the gestational length in high-
risk pregnancies. Fundamentally, these questions can be understand as problems
of binary classification (is a pregnancy at risk of being preterm?), anomaly
detection (can we identify high risk pregnancy when the prevalence is “only”
5-18%?), and survival analysis (can we understand gestational length as a
time-to event?). The ability to answer these questions promises to improve the
personalization of prenatal care.

Recently, the microbiome has emerged in the study of preterm delivery and one
dataset that is uniquely suited to addressing all three questions of
classification, anomaly detection, and survival analysis is from DiGiulio et
al. (2015). The authors tracked the microbiomes of 40 women over the course of
their pregnancies by collecting over 3400 samples from the vagina, distal gut,
saliva, and tooth/gum areas. This effort resulted in a rich dataset that
revealed the presence of over 1200 different taxa. The authors’ aim was to
identify microbial taxa associated with a higher risk of preterm delivery,
which occurred in 15 of the 40 women in the study. We summarize the lengths of
gestation in Figure 1.

Figure 1: Histogram of lengths of gestational periods in weeks. The dashed
line separates the classes of term, which is delivery after 37 weeks, from
preterm. Red and blue represent women with and without a history of preterm
delivery, respectively.

In their analysis, DiGiulio et al. (2015) successfully identified several taxa
whose dynamics throughout the gestational period were significantly related to
a higher risk of preterm delivery. Additionally, these authors looked at the
interconnectedness of taxa at the level of taxonomic communities. Ultimately,
corresponding network summary values were used as descriptors for linear
mixed-effect models. The resulting findings suggest that individual
microbiomes – as captured through network-based representations – can be
important descriptors for various pregnancy outcomes. For a recent review of
the current state of applying network analyses to microbiome data, see
Layeghifard, Hwang and Guttman (2017).

Examples of such networks are shown in Figure 2. In DiGiulio et al. (2015) and
other microbiome analyses, the strategy has been to identify and extract
certain potentially relevant aspects of the topology of these networks for
downstream analysis. However, it is clear from these visualizations that
microbiome networks possess a rich topology and that they can exhibit
important topological differences across, for example, patient subgroups.
Without knowing precisely which aspects of network topology might be most
relevant to a given study, it is increasingly desirable to have statistical
methods that allow researchers to incorporate the entire network, rather than
just numerical summaries.

Figure 2: Examples of individual microbiome networks. Nodes represent
different species of bacteria. They are colored by phylum and their size is
proportional to their abundance across the individual’s samples. Edges
correspond to association between the species with the edge thickness
corresponding to the strength of association and color corresponding to the
sign (black is negative and red is positive). Further details of the network
constructions are outlined in Section 6.1. Although it appears that the
preterm microbiome network has a lower density and is dominated by the
presence of a few species, we show in Section 6.2 that the differences in
preterm and term microbiome networks cannot be easily reduced to a few summary
statistics.

Herein, we propose a unified methodology for analyzing network datasets for a
number of standard purposes using whole-network inputs. In particular, we
develop a classifier that simultaneously serves as the basis for the
corresponding anomaly detection and survival analysis methodologies we also
develop in this paper in a seamless and integrated manner. A schematic diagram
of the core infrastructure is given in Figure 3.

We expect our contributions to advance the use of networks to analyze the
microbiome, both now and in the future. Ultimately, a better understanding of
the relationship between preterm delivery and microbiome networks that are
sufficiently easy to obtain and construct could help to inform routine
prenatal care decisions, as well as to provide important controls in studies
of pregnant women. Specifically, our Bayesian approach can help mitigate the
challenges of small sample size (McNeish, 2016), which typifies most current
microbiome experiments, and at the same time our approach scales to the
expected paradigm shift toward microbiome studies that sample with higher
resolution for many individuals.

Figure 3: Schematic diagram of study design, network construction, and
probabilistic modeling with network inputs. The $4\times 4$ symmetric matrix
$K$ represents a kernel that is used to measure (dis)similarity between pairs
of networks. In this toy example, the color entries are pairwise averages of
the colors of the corresponding networks, e.g. the purple in entry $K_{1,2}$
is the average of the blue network (row 1) and the orange network (column 2).

###  1.1 Background

We briefly review relevant network and Gaussian process (GP) concepts. For
more thorough reviews of the statistical analysis of network data and the use
of GPs in machine learning, see Kolaczyk and Csárdi (2014) and Rasmussen and
Williams (2005), respectively.

A network or graph $G$ is an ordered pair $(V,E)$ consisting of a set of nodes
or vertices, $V$, and edges, $E$, where $e\in E(G)$ is a pair of nodes
$(v_{1},v_{2})\in V(G)\times V(G)$. Throughout, we work with undirected and
simple networks, meaning that there is no directional component to our edges
and no self loops. Furthermore, we will work with networks of the same order,
i.e. for networks $G_{1},\ldots,G_{m}$, we have
$|V_{1}|=|V_{2}|=\cdots=|V_{m}|=n$.

For computations and analyses, graphs are represented through matrices. The
adjacency matrix of a graph $G$ with nodes $V(G)=(v_{1},\ldots,v_{n})$ is an
$n\times n$ matrix $A$ with

$A_{ij}=\begin{cases}1&\text{if }v_{i}\text{ adjacent to }v_{j}\\\
0&\text{otherwise}\enskip.\end{cases}$

For weighted networks, which are triples $G=(V,E,W)$, the corresponding
weighted adjacency matrix is defined similarly:

$A_{ij}=\begin{cases}w_{ij}&\text{if }\\{i,j\\}\in E\\\
0&\text{otherwise}\enskip.\end{cases}$

A Gaussian process (GP) is a stochastic process that is commonly used for
modeling a random function whose evaluation at a finite number of points is
jointly Gaussian. In particular, a GP is often used as a prior distribution
for the underlying regression function or latent classification function in a
Bayesian nonparametric model which leads to a posterior distribution of the
function for inference. We write

$f(x)\sim\mathcal{GP}\big{(}m(x),\ k(x,x^{\prime})\big{)}\enskip,$

where $m(x)=\mathbb{E}[f(x)]$ is the mean function and
$k(x,x^{\prime})=\mathbb{E}\big{[}\big{(}f(x)-m(x)\big{)}\big{(}f(x^{\prime})-m(x^{\prime})\big{)}\big{]}$
is the covariance kernel. Then for any $x_{1},\ldots,x_{n}$,
$\left(f(x_{1}),\ldots,f(x_{n})\right)$ is multivariate Gaussian with mean
vector $\left(m(x_{1}),\ldots,m(x_{n})\right)$ and covariance matrix
$K=((k(x_{i},x_{j}))_{1\leq i\leq j\leq n}$. One of the most widely used
covariance kernels for $x\in\mathbb{R}^{d}$ is the squared-exponential kernel,

$k(x_{i},x_{j})=\sigma^{2}\exp\Big{(}-\frac{|x_{i}-x_{j}|^{2}}{2\ell^{2}}\Big{)}\enskip,$
(1)

where $|\cdot|$ refers to Euclidean distance and $\sigma^{2}$ and $\ell$ are
the signal variance and length-scale parameters.

Note that GPs can be defined on any topological space. In our case, this
corresponds to the space of all networks on $n$ nodes. Therefore, each $x$
corresponds a network input. We propose to use a GP as a prior distribution
for modeling a function on networks which could correspond to a regression
function, latent classification map, or the survival curve in a given
statistical model with networks given as covariate information. One of the key
difficulties in using a GP model for networks lies in constructing a valid
covariance kernel which is a positive definite kernel on the space of
networks. We propose several kernels in Section 2.1 that aim to capture
different networks properties.

The covariance kernel such as the squared-exponential type given in (1) is
equipped with additional hyperparameters that need to be learned. To do so, we
adopt a Bayesian framework by employing hierarchical models for our parameters
and hyperparameters. This has the further benefit of being probabilistic,
which, in addition to providing uncertainty quantification of our posterior
predictions, yields an immediate extension from our classifier to anomaly
detection. Finally, the Bayesian framework allows for prior specification such
as clinicians informing priors for the baseline hazard function in a survival
analysis.

###  1.2 Related Work

There is a large literature on graph kernels beginning with Kondor and
Lafferty (2002), who propose diffusion kernels for graphs and suggest that
they could be used in conjunction with GPs. Shortly thereafter, Gärtner,
Driessens and Ramon (2003) propose the random walk graph kernel, which the
authors use to perform a variant of GP regression. To the best of our
knowledge, this is the only explicit use of GPs with graphs. Kashima and
Inokuchi (2002) apply similar graph kernels for classification, but not using
a GP framework. Since then, there have been many other graph kernels proposed.
For surveys on graph kernels, see Vishwanathan et al. (2010); Kriege,
Johansson and Morris (2020); Nikolentzos, Siglidis and Vazirgiannis (2019).

Regarding the use of network inputs, a Bayesian approach seems not to have
been applied to classification problems. However, kernel support vector
machines (SVMs) have been a popular tool for classification with network
inputs and extensions exist to one-class classification and survival analysis.
In particular, graph data has been used with SVMs to perform classification of
protein function prediction (Borgwardt et al., 2005), chemical informatics
(Ralaivola et al., 2005), and disease (Rudd, 2018), as well as one-class
classification for media data (Mygdalis et al., 2016). For the connection
between the objective functions of SVMs and GP classification, see (Rasmussen
and Williams, 2005, Chapter 6).

Alternatively, Arroyo Relión et al. (2019) provide a frequentist approach to
classification with network inputs. Interestingly, in comparing their method
to others, the authors say that “kernel methods were no better than random
guessing.” They also claim that “kernel methods [are] unsuitable for large
scale networks.” Throughout, we demonstrate that, in fact, our method is
highly scalable to large networks and outperforms those authors’ method on
their application and other classification tasks.

Elsewhere, Zhang et al. (2013) introduce Net-Cox for survival analysis with
network inputs, which is a network-constrained Cox regression using the graph
Laplacian as a penalty. The authors note that “surprisingly, network-based
survival analysis has not received enough attention.” We will discuss the
benefits of our survival analysis methodology, which is not a proportional
hazards model, in Section 2.3.

Finally, there is a growing body of work on statistical analysis of multiple
networks in general. For example, Jain and colleagues have a number of
contributions, summarized, for example, in Jain (2016), where they propose the
use of linear classifiers in graph space. While our work also has some
geometrical underpinnings, it is nevertheless distinct. Similarly, there is
various work on network-based extensions of averages (Ginestet et al., 2017;
Durante, Dunson and Vogelstein, 2017; Tang et al., 2018; Kolaczyk et al.,
2020), regression (Cornea et al., 2017), clustering (Mukherjee, Sarkar and
Lin, 2017), and PCA (Dai et al., 2018), to name a few. Finally, there are
several latent space models for multiple networks, for example in Arroyo et
al. (2019); Gollini and Murphy (2016); Salter-Townshend and McCormick (2017).
Again, however, our work is distinct.

###  1.3 Paper Outline

The remainder of the paper is organized as follows. In Section 2, we set up
our models for binary classification, anomaly detection, and survival
analysis. We discuss implementation and theoretical results for our models in
Section 3 and Section 4, respectively. We present the results of our models on
simulated networks in Section 5, as well our analysis of the microbiome data
in Section 6. Finally, we conclude in Section 7 with a discussion of future
directions for this work.

##  2 Models

###  2.1 Binary Classification

Suppose we have data $(G_{1},Y_{1}),\ldots,(G_{m},Y_{m})$, where

$G\in\mathcal{G}:=\big{\\{}\text{(labeled) weighted networks of order
}n\big{\\}}$

and $Y$ is binary. We are interested in learning the classification map
$\pi(G):=p(y=1\ |\ G)$. To do this, we will directly model $\pi(G)$ using
Bayesian inference. In particular, we equip $\pi(G)$ with a prior that is a
deterministic transformation of a GP on $\mathcal{G}$. That is,
$\pi(G)=H(f(G))$ for some link function $H$ and latent function $f:=f(G)$.
Throughout, we work with the logistic link,
$\text{logistic}(z)=1/(1+\exp(-z))$. We then place a GP prior on $\mathcal{G}$
with kernel $K$ over $f$, where $K$ is given by

$K_{ij}=\sigma^{2}\exp\Big{(}-\ell\cdot d(G_{i},G_{j})\Big{)}\enskip.$ (2)

The squared-exponential kernel in (2) is infinitely differentiable and
stationary, which we will use to show posterior consistency in Theorems 4.2
and 4.3. The challenge in working with this kernel is in specifying an
appropriate graph distance. Although there are many such distances (Donnat et
al., 2018), a squared-exponential kernel is positive definite if and only if
its distance is induced by an inner product (Jayasumana et al., 2013), in this
case on the space of networks $\mathcal{G}$. In Section 4.1, we show that the
following distances satisfy this requirement:

$d_{\text{F}}\big{(}G_{1},G_{2}\big{)}:=\frac{1}{n(n-1)}\sum_{i,j\in
V}\big{(}A_{ij}^{(1)}-A_{ij}^{(2)}\big{)}^{2}\enskip,$ (3)

and

$d_{\lambda}\big{(}G_{1},G_{2}\big{)}:=\sum_{i=0}^{n-1}\big{|}\lambda_{i}^{(1)}-\lambda_{i}^{(2)}\big{|}^{2}\enskip,$
(4)

where $A^{(k)}$ is the adjacency matrix for network $G_{k}$ and
$\lambda_{0}^{(k)}\leq\lambda_{1}^{(k)}\leq\cdots\leq\lambda_{n-1}^{(k)}$ are
the eigenvalues of its Laplacian.

The distance in (3) is the Frobenius distance, a slight modification of the
Hamming distance for binary networks, which captures local measure of
difference, and the distance in (4) is a spectral-based distance with a more
global measure of difference.

We note that graph-based distances may not always capture the complex topology
of the network’s distribution. There are many other graph kernels proposed in
the literature, but the fundamental challenge is ensuring its positive
definiteness. One positive definite alternative is the random-walk kernel in
Vishwanathan et al. (2010), which can be easily plugged in to our method,
though it has fewer kernel hyperparameters. In Section 5, we isolate the
benefits of our GP approach by comparing a GP classifier with a random walk
kernel against an SVM with the same kernel.

For a squared-exponential kernel, our model, including inverse-gamma priors on
the kernel hyperparameters and the posterior from which we will make inference
on $\pi(G)$, is specified as follows:

likelihood $\displaystyle:\hskip 7.11317ptp(y\ |\
f)=\prod_{k=1}^{m}\text{logistic}(f_{k}y_{k})\enskip;$ prior
$\displaystyle:\hskip 7.11317ptp(f\ |\ G,\sigma^{2},\ell)\propto
e^{-\frac{1}{2}f^{T}K^{-1}f}\enskip,$ $\displaystyle\ \ \ \hskip
7.11317ptp(\sigma^{2})\propto(\sigma^{2})^{-\alpha_{1}-1}e^{-\frac{\beta_{1}}{\sigma^{2}}}\enskip,$
$\displaystyle\ \ \ \hskip
7.11317ptp(\ell)\propto\ell^{-\alpha_{2}-1}e^{-\frac{\beta_{2}}{\ell}}\enskip;$
posterior $\displaystyle:\hskip 7.11317ptp(f,\sigma^{2},\ell\ |\ G,y)\propto
p(y\ |\ f)p(f\ |\ G,\sigma^{2},\ell)p(\sigma^{2})p(\ell)\enskip.$

Again, any positive definite kernel may be used in place of $K$ and the only
difference would be the prior specification of its hyperparameters (or removal
of these priors, for example in the case of the random walk kernel).

Regarding our likelihood, note that if we take $Y\in\\{\pm 1\\}$, then we can
write

$p(y\ |\ f)=\text{logistic}(f\cdot y)\enskip,$

as $\text{logistic}(-z)=1-\text{logistic}(z)$ and $p(y=1\ |\ G)+p(y=-1\ |\
G)=1$.

For a new network $\tilde{G}$, the posterior predictive distributions are

$\displaystyle p(\tilde{f}\ |\ G,y,\tilde{G})$ $\displaystyle\propto
p(\tilde{f}\ |\ G,\tilde{G},f)p(f\ |\ G,y)\enskip,$ $\displaystyle
p(\tilde{y}=1\ |\ \tilde{G},G,y)$
$\displaystyle=\int\text{logistic}(\tilde{f})\cdot p(\tilde{f}\ |\
G,y,\tilde{G})d\tilde{f}\enskip,$

where

$p(\tilde{f}\ |\ G,\tilde{G},f)\sim
N\Big{(}k(\tilde{G})^{T}K^{-1}f,k(\tilde{G})-k(\tilde{G})^{T}K^{-1}k(\tilde{G})\Big{)}\enskip.$

Not surprisingly, our posterior is intractable, which means that our posterior
predictive distributions are also intractable. In Section 3.1, we discuss
implementation strategies including Markov chain Monte Carlo (MCMC) sampling
and Monte Carlo integration.

###  2.2 Anomaly Detection

One-class classification (OCC) is a classifier for a single class. For a
review, see Khan and Madden (2009). Such a classifier is a type of anomaly
detection that is useful when most of the data is from a single class and only
a few, if any, training points are from another class.

The most common approach to OCC is through decision boundaries, which are
typically found via SVMs. There are many variations of one-class SVM (OSVM)
depending on the availability of training data, i.e whether any training data
is from the negative class. Another approach to OCC is to simply use machine
learning techniques for classification, like k-NN or tree-based methods, but
these require negative training examples. However, neither OSVM nor non-OSVM
solutions are probabilistic.

GP methods have been proposed as well for the OCC problem. Kemmler et al.
(2013) suggests four GP values for OCC scores:

$\displaystyle\tilde{\mu}$ $\displaystyle=\mathbb{E}[\tilde{f}\ |\
G,y,\tilde{G}]\enskip,$ $\displaystyle\tilde{\pi}$
$\displaystyle=p[\tilde{y}=1\ |\
G,y,\tilde{G}]=\text{logistc}(\tilde{\mu})\enskip,$
$\displaystyle\tilde{\sigma}^{2}$ $\displaystyle=\text{Var}[\tilde{f}\ |\
G,y,\tilde{G}]\enskip,$ H
$\displaystyle=\tilde{\mu}\cdot\tilde{\sigma}^{-1}\enskip.$

These values are obtained immediately from the posterior of our classification
model. Unsurprisingly, Kemmler et al. (2013) notes that “tuning
hyperparameters for one-classification tasks is a difficult task in a general
setting without incorporating further model assumptions,” but we overcome this
challenge through our Bayesian approach.

###  2.3 Survival Analysis

Fernández, Rivera and Teh (2016) introduce a semi-parametric Bayesian method
for survival analysis that uses GPs to model variation around a parametric
baseline hazard function. Their method easily incorporates covariates,
censoring, and prior knowledge, while avoiding the proportional hazards
constraint. By using our kernel inputs, we are able to adapt this model to
perform survival analysis with network inputs.

Suppose we have data $(G_{1},T_{1}),\ldots,(G_{m},T_{m})$, where
$G\in\mathcal{G}$ and $T$ is a survival time on $\mathbb{R}^{+}$ with survival
function $S$ and hazard function $\lambda$. The authors model $T$ as the first
jump of a Poisson process with intensity $\omega$. Together with the GP
priors, the model for $T_{i}\ |\ G_{i}$ is

$\displaystyle T_{i}\ |\ \omega_{i}$
$\displaystyle\overset{ind}{\sim}\omega(T_{i})e^{-\int_{0}^{T_{i}}\omega_{i}(s)ds}\enskip,$
$\displaystyle\omega_{i}(t)\ |\ f,\omega_{0}(t),G_{i}$
$\displaystyle\overset{\hphantom{ind}}{=}\omega_{0}(t)\cdot\text{logistic}(f(t,G_{i}))\enskip,$
$\displaystyle f(\cdot)$
$\displaystyle\overset{\hphantom{ind}}{\sim}\mathcal{GP}(0,K)\enskip,$

where $K$ is a kernel in network and time. We take

$K\big{(}(t_{1},G_{1}),(t_{2},G_{2})\big{)}=K_{T}(t_{1},t_{2})+K_{G}(G_{1},G_{2})\enskip,$
(5)

where $K_{T}$ and $K_{G}$ are the kernels given in is a kernel for time given
in (1) and (2), respectively, with shared signal variance $\sigma^{2}$. The
authors prove (Proposition 1) that for stationary kernels such as the squared-
exponential, $S(t)$ associated with $f(t)$ is a proper survival function, i.e.
for a fixed network $G$, we have $S_{G}(t)=\mathbb{P}(T>t\ |\ G)\to 0$ as
$t\to\infty$.

Unfortunately, the likelihood for $T_{i}\ |\ G_{i}$ is (doubly) intractable
since $\omega_{i}$ is defined by a GP. To overcome this, the authors develop a
data augmentation scheme for sampling from an inhomogeneous Poisson process
with intensity $\omega_{0}(t)$, which allows a tractable reformulation of the
model. Using the tractable model, which is a joint distribution on $(R,T)$,
where $R$ are (unknown) rejected jump points from a thinned Poisson process
and $T$ is the (known) first accepted one, the authors develop an inference
algorithm that begins by sampling $R$. Crucially, the inference algorithm
relies on sampling the GP $f$ given $R\cup T$, which can be understood as GP
binary classiﬁcation. We provide implementation details in Section 3.2,
including remarks on how our implementation differs from the original based on
our use of network inputs.

##  3 Implementation

In this section, we provide details for implementing the models from Section
2. All of the  code for our algorithms – along with the code for reproducing
our figures, simulations, and data analysis – is available at
github.com/KolaczykResearch/GP-Networks.

###  3.1 Binary Classification and OCC

We develop a Gibbs sampler to overcome the intractability of our posterior.
The posterior conditional distributions are as follows:

$\displaystyle\sigma^{2}\ |\ f,\ell$ $\displaystyle\sim\text{Inv-
Gamma}\Big{(}\alpha_{\sigma}+\frac{m}{2},\
\beta_{\sigma}+\frac{f^{T}K_{0}^{-1}f}{2}\Big{)}\enskip,$ $\displaystyle\hskip
7.11317ptp(\ell\ |\ f,\sigma^{2})$
$\displaystyle\propto|K|^{-1/2}e^{-\frac{1}{2}f^{T}K^{-1}f-\frac{\beta_{\ell}}{\ell}}\ell^{-\alpha_{\ell}-1}\enskip,$
$\displaystyle p(f\ |\ \sigma^{2},\ell)$ $\displaystyle\propto
e^{-\frac{1}{2}f^{T}K^{-1}f}\prod_{k=1}^{m}\text{logistic}(f_{k}y_{k})\enskip.$

The signal variance parameter is conjugate, which makes sampling easy and in
general this hyperparameter just controls the scaling of our latent function.
However, the length scale parameters are coupled with the latent function,
which makes sampling more challenging. Moreover, the determinant term adds
unnecessary computational time if we elect to sample $f$ and $\ell$
separately. Instead, we sample them jointly following Murray and Adams (2010)
Additionally, we make use of an elliptical slice sampler (Murray, Adams and
MacKay, 2010) for better mixing of the latent function, which adds essentially
no cost as the Cholesky is already cached. The pseudocode for our Gibbs
sampler is given in Algorithm 1.

GPC Gibbs $(K,Y,\alpha,\beta,ns)$;

Input : Kernel $K$, Labels $Y\in\\{\pm 1\\}$, inverse-gamma parameters
$\alpha\ \&\ \beta$, Number of samples $ns$

Output : Posterior $p(f|y,\theta)$

Initialize $\sigma^{2},\ell,f$;

for  _$t=2,\ldots,ns$_ do

Jointly sample $\big{(}f^{(t)},\ell^{(t)}\big{)}\ |\ \sigma^{(t-1)}$ using
slice sampler;

Compute $C_{0}$, the Cholesky of $K_{0}$ evaluated at $\ell^{(t)}$;

$C=C_{0}\times\sqrt{\sigma^{(t-1)}}$, the Cholesky of $K$ evaluated at
$\sigma^{(t-1)},\ell^{(t)}$;

Sample $f^{(t)}\ |\ \sigma^{(t-1)},\ell^{(t)}$ cheaply using elliptical slice
sampler;

Sample $\sigma^{2}_{t}\sim\text{Inv-
Gamma}\big{(}\alpha_{\sigma}+\frac{m}{2},\beta_{\sigma}+\frac{f^{T}K_{0}^{-1}f}{2}\big{)}$;

$C=C_{0}\times\sqrt{\sigma^{(t)}}$

end for

Algorithm 1 Gibbs algorithm for GP classification.

Note that the complexity of our sampler is dominated by the Cholesky
decomposition of our kernel, which is a limitation inherent to GP methods.
However, contrary to the claim in Arroyo Relión et al. (2019), our method
scales well in the size of the network, because a distance matrix $D$, with
$D_{ij}=d(G_{i},G_{j})$, only needs to be computed once and can be done in
parallel. Of course, other implementation strategies such as a grid search or
maximizing the marginal likelihood would be faster than ours. However, these
approaches lose the full benefit of a Bayesian framework and in practice, we
find they perform worse, as they are more limited in learning the kernel
hyperparameters.

Suppose we have run Gibbs on a training set. To make predictions, first draw
$B$ samples from the posterior. Then, we can estimate $\tilde{f}$ in two ways:

$\displaystyle\hat{\tilde{f}}$
$\displaystyle=\frac{1}{B}\sum_{b=1}^{B}k_{\theta^{(b)}}(\tilde{G})^{T}K_{\theta^{(b)}}^{-1}\hat{f}\enskip,$
(6) $\displaystyle\hat{\tilde{f}}_{\text{avg}}$
$\displaystyle=k_{\hat{\theta}}(\tilde{G})^{T}K_{\hat{\theta}}^{-1}\hat{f}\enskip,$
(7)

where $\theta^{(b)}$ is the $b^{\text{th}}$ sample of
$\theta=(\sigma^{2},\ell)$ and $\hat{f}$ and $\hat{\theta}$ are the estimated
posterior means of $f$ and $\theta$, respectively, all of which come from the
Gibbs sampler. We recommend the plug-in estimate via the posterior mean in
(7), which only requires one matrix inversion, as in practice we see little
difference in accuracy between the versions.

###  3.2 Survival Analysis

The posterior conditional distributions are as follows:

$\displaystyle\sigma^{2}\ |\ f,\ell_{T},\ell_{\text{G}}$
$\displaystyle\sim\text{Inv-Gamma}\Big{(}\alpha_{\sigma}+\frac{m+|R|}{2},\
\beta_{\sigma}+\frac{f^{T}K_{0}^{-1}f}{2}\Big{)}\enskip,$ $\displaystyle
p(\ell_{p}\ |\ f,\sigma^{2})$
$\displaystyle\propto|K|^{-1/2}e^{-\frac{1}{2}f^{T}K^{-1}f-\frac{\beta_{p}}{\ell_{p}}}\ell_{p}^{-\alpha_{p}-1}\enskip,$
$\displaystyle p(f\ |\ \sigma^{2},\ell)$ $\displaystyle\propto
e^{-\frac{1}{2}f^{T}K^{-1}f}\prod_{k=1}^{m}\text{logistic}\big{(}f(T_{k})\big{)}\prod_{r\in
R_{k}}\big{(}1-\text{logistic}(f(r)\big{)}\enskip.$

The pseudocode is given in Algorithm 2.

GP Survival $(D,T,\alpha,\beta,ns)$;

Input : Distance array $D$, Survival times $T$, Hyperparameter priors $\alpha\
\&\ \beta$, Number of samples $ns$

Output : Samples from $\omega\ |\ T$

Initialize baseline hazard $\omega_{0}$ and kernel $K$;

$\lambda_{0}=\Omega$ and $\Lambda_{0}(t)=\Omega\cdot t$;

$K_{ij}=\sigma^{2}\Big{(}\exp(-\ell_{G}\cdot
D_{ij}^{G})+\sum_{k=1}^{p}\exp(-\ell_{k}\cdot D_{ij}^{p})+\exp(-\ell_{T}\cdot
D_{ij}^{T}\Big{)}$;

Instantiate $f$ in $T$;

$f^{(1)}(T)\sim\mathcal{N}(0,K)$;

for  _$t=2,\ldots,ns$_ do

for  _$i=1,\ldots,m$_ do in parallel

$n_{i}\sim\text{Poisson}(1;\Lambda_{0}(T_{i}))$;

$\tilde{A}_{i}\sim\mathcal{U}(n_{i};0,\Lambda_{0}(T_{i}))$;

$A_{i}=\Lambda_{0}^{-1}(\tilde{A}_{i})$;

Sample $f(A_{i})\ |\ f(R\cup T),\lambda_{0}$;

$f^{(t)}(A_{i})\sim\mathcal{N}\big{(}k(A_{i},R\cup T)^{T}K^{-1}f^{(q-1)}(R\cup
T),\ k(A_{i},A_{i})-k(A_{i},R\cup T)K^{-1}K(R\cup T,A_{i})\big{)}$;

$U_{i}\sim\mathcal{U}(n_{i};\ 0,1)$;

$R_{i}=\\{a\in A_{i}\text{ such that }U_{i}<1-\sigma(f(a))\\}$;

end forpar

$R=\bigcup_{i=1}^{m}R_{i}$;

Update parameters of $\lambda_{0}$;

$\Omega\sim\Gamma\big{(}\alpha_{\Omega}+m+|R|,\beta_{\Omega}+\sum_{i=1}^{m}T_{i}\big{)}$;

Update $f(R\cup T)$ and hyperparameters of kernel as in Algorithm 1;

Jointly sample $f(R\cup T)$ and length scales for time and covariates using
slice sampler;

Sample $f(R\cup T)$ several times using elliptical slice sampler;

Sample signal variance from conjugate posterior;

end for

Algorithm 2 Gibbs algorithm for GP survival analysis with constant baseline
hazard function.

Note that $f$ changes size at each iteration of the sampler: at $t=1$,
$R=\emptyset$ and $|f|=m$, whereas for $t\geq 2$, we have $|f|=m+|R|$. To
update $f(R\cup T)$ in our slice sampler, we input the concatenation of $f(T)$
with $f(R)$, where the latter is from our previous sample $f(A)$. Also, notice
that the sampling of rejected points and their imputation given $f(R\cup T)$
can be implemented in parallel, which was not discussed in the original paper
(Fernández, Rivera and Teh, 2016). In addition, it is important to remember
that each time point has associated covariates that are omitted in the
notation, i.e. $f(t)=f(t,X)$, where $X$ in this case is a network. Therefore,
it is necessary to keep track of the implicit covariates when sampling $A$.
Finally, as the authors point out, setting $R$ can be seen as a GP
classification problem. This procedures essentially uses a noisy version of
our classifier from Section 2.1, which would be too computationally expensive
to use at every step of the Markov chain.

One computational note is that Algorithm 2 can be slow if $|R|$ is large,
which the authors circumvent using a kernel approximation. Unfortunately, such
approximations are not available for network kernels because they rely on
Bochner’s theorem and Fourier transforms of kernels on $\mathbb{R}^{d}$.
Another computational note is that the kernels, while theoretically positive
definite, can be numerically close to singular. For this reason, it is common
to add jitter to the diagonal of the kernel to improve the condition number.
We find that more jitter is needed in the survival analysis setting than the
binary classification, which is likely due to the fact that the rejected
points $R$ share all of the covariates with one of the original points,
thereby introducing dependency into the rows and columns of the Gram matrix.

Having run the sampler, the survival surfaces are given as

$S(t,X)=\exp\Big{(}-\int_{0}^{t}\omega(s,X)\ ds\Big{)}\enskip.$

We use the trapezoidal rule for integration by evaluating $f(s,X)$ on a
uniform grid from $s=0$ to $s=\max T$ with $\Delta=\frac{\max T}{99}$, i.e.

$\int_{0}^{t}\omega(s,X)\
ds\approx\frac{\Delta}{2}\sum_{k=1}^{100}\omega(s_{k-1},X)+\omega(s_{k},X)\enskip.$

##  4 Theoretical Results

In this section, we prove that our graph distances in (3) and (4) induce valid
kernels. Using these kernels, we then prove that our binary classifier from
Section 2.1 and our survival analysis model from Section 2.3 both achieve
posterior consistency.

###  4.1 PD Kernels

Central to the GP framework is a kernel for capturing the dissimilarity
between data points. However, not all dissimilarity measures will induce a
valid kernel, i.e a provably positive definite Gram matrix. This problem is
further complicated by our use of network objects as it is not obvious which
graph distances will ensure a valid kernel.

In order to work with the squared-exponential kernel, we need a distance that
is an inner product on the space of graphs. One familiar choice is the Hamming
distance or Frobenius distance, which is typically defined for two binary
networks as the sum of the absolute differences of their corresponding
adjacency matrices. However, as we will see in our proof of consistency, we
need to work in the space of weighted networks. Moreover, we often encounter
weighted networks, so using the edge weights rather than indicators for edges
may provide us more information and we avoid the common challenge of
determining a threshold value for edges. In this setting, we can view binary
networks as having weights $\\{0,1\\}\subset[0,1]$.

In order to work with weighted networks, we use the Frobenius distance in (3),
which we recall as

$d_{\text{F}}\big{(}G_{1},G_{2}\big{)}:=\frac{1}{n(n-1)}\sum_{i,j\in
V}\big{(}A_{ij}^{(1)}-A_{ij}^{(2)}\big{)}^{2}\enskip.$

Note that for binary networks, the sum of absolute differences is equal to the
sum of squared differences, so our Frobenius distance is consistent with the
commonly used Hamming distance for binary networks. Furthermore, the Frobenius
distance can be used for graphs with negative weights such as measures of
correlation.

Another choice of graph distance is in (4), which we recall as

$d_{\lambda}\big{(}G_{1},G_{2}\big{)}:=\sum_{i=0}^{n-1}\big{|}\lambda_{i}^{(1)}-\lambda_{i}^{(2)}\big{|}^{2}\enskip.$

This is a spectral-based distance, which has the benefit of capturing more of
the network structure at the cost of more computational complexity.

######  Theorem 4.1.

If $G_{1},\ldots,G_{m}$ are weighted networks and

$K_{ij}=\sigma^{2}\exp\Big{(}-\ell\cdot
d\big{(}G_{i},G_{j}\big{)}\Big{)}\enskip,$

with $d(\cdot,\cdot)$ given in (3) or (4), then $K$ is a positive-definite
kernel.

###### Proof.

Jayasumana et al. (2013) provide the following theorem. Let $(M,d)$ be a
metric space and define $k:X\times X\to\mathbb{R}$ as
$k(x_{i},x_{j})=\exp\big{\\{}-d^{2}(x_{i},x_{j})/2\sigma^{2}\big{\\}}$. Then
$k$ is positive definite for all $\sigma>0$ if and only if there exists an
inner product space $V$ and a function $\psi:M\to V$ such that
$d(x_{i},x_{j})=\|\psi(x_{i})-\psi(x_{j})\|$.

Using the Frobenius distance and absorbing the normalizing constant into
$\sigma^{2}$, we have

$\displaystyle d_{\text{F}}\big{(}G_{1},G_{2}\big{)}$
$\displaystyle=\sum_{i,j\in V}\big{(}A_{ij}^{(1)}-A_{ij}^{(2)}\big{)}^{2}$
$\displaystyle=(A^{(1)}-A^{(2)})\cdot(A^{(1)}-A^{(2)})$
$\displaystyle=\|A^{(1)}-A^{(2)}\|_{\text{F}}^{2}\enskip,$

where $\cdot$ and $\|\cdot\|_{\text{F}}$ are the Frobenius inner product and
norm, respectively, and $\psi(\cdot)$ is the identity map.

Similarly, using the spectral-based distance, we have

$d_{\lambda}\big{(}G_{1},G_{2}\big{)}=\sum_{i=0}^{n-1}\big{|}\lambda_{i}^{(1)}-\lambda_{i}^{(2)}\big{|}^{2}=\|\lambda^{(1)}-\lambda^{(2)}\|^{2}\enskip,$

where $\psi(\cdot)$ is the map taking the graph’s adjacency matrix to the
eigenvalues of its Laplacian. ∎

###  4.2 Posterior Consistency of Classifier

Let $G_{1},\ldots,G_{m}$ be simple weighted networks on the same set of $n$
vertices, $V$. That is, $G_{k}=(V,E^{(k)},W^{(k)})$ for all $k=1,\ldots,m$
with $w_{jj}=0$ for all $j=1,\ldots,n$. Unlike many graph results that
consider the asymptotic behavior as the order of the network, $n$, grows, we
return to the classical data setting and investigate what happens as the
sample size, $m$, of observed networks increases. We expect that this
asymptotic regime will become increasingly important in the network literature
as it is likely, for example in the case of connectomics, that we sample more
brain networks with a fixed number of regions of interest.

With the set up from Section 2.1 with inverse-gamma hyperpriors over the
kernel induced by either distance (3) or (4), we show that our classifier
achieves posterior consistency.

######  Theorem 4.2.

GP classification with network inputs achieves posterior consistency for the
squared-exponential kernel.

We defer the proof to the Appendix. Ultimately, we show that our classifier
achieves posterior consistency by verifying the conditions given in Ghosal and
Roy (2006). Here, we make a few remarks about the proof.

First, the infinite differentiability of our kernel is important for verifying
several of the conditions, which demonstrates the importance of using a
squared-exponential kernel rather than, say, a Matérn kernel. This reiterates
the need to have a graph distance that guarantees positive-definiteness of our
squared-exponential kernel. Secondly, it is important that we bound our
covariate space, which forces us to work in the space of weighted networks.
Finally, the inverse-gamma hyperpriors ensure full support over the kernel
hyperparameters. Although other fully supported priors could be used while
maintaining posterior consistency, inverse-gamma is conjugate for the signal
variance hyperparameter.

###  4.3 Posterior Consistency of Survival Analysis

As above, let $G_{1},\ldots,G_{m}$ be simple weighted networks on the same set
of vertices. Consider the survival analysis model from Section 2.3 with the
kernel in network and time given by (5). Under the assumption that the
survival times have finite expectation, we show that our survival analysis
model achieves posterior consistency in the limit of $m$.

######  Theorem 4.3.

Under mild assumptions, survival analysis with network inputs achieves
posterior consistency.

Again, we defer the proof to the Appendix, which consists of verifying the
conditions given in Fernández and Teh (2016). The main assumption that we make
is that the mean of the survival time is finite, which is true for events like
death or gestational age and reasonable for many others. Otherwise, similar to
Theorem 4.3, we rely on our use of a squared-exponential kernel due to its
stationarity.

##  5 Simulations

In this section, we evaluate our models for binary classification, anomaly
detection, and survival analysis in several simulated settings.

###  5.1 Binary Classification

We first present the results of a simulation study to assess the accuracy of
our GP classifier. Since classification is at the center of our anomaly
detection and survival analysis models, it is important to identify what
scenarios are better suited to one or another kernel or distance.

We consider the squared-exponential kernel with both the Frobenius distance in
(3) and the spectral-based distance in (4) as well as a $k$-step random walk
kernel in Vishwanathan et al. (2010). We refer to these methods as GP-F,
GP-$\lambda$, and GP-RW, respectively. To isolate the contribution of our
method, we also include a comparison to an SVM with the same random walk
kernel. We use a cross-validation procedure to choose the cost penalty of the
SVM from the grid suggested by Hsu et al. (2003) and we refer to this method
as SVM-RW. Finally, we also include two recent graph classifiers from Arroyo
Relión et al. (2019) and Arroyo et al. (2019), which we refer to as GC and
MASE.

Our simulation is a fully crossed design varying the sample size
($m=40,100,\ldots,340$) and the network size ($n=20,40,\ldots,100$), as well
as the random graph model from which we sample the network inputs. We compare
small-world networks (Watts and Strogatz, 1998) with different rewiring
probabilities $p0=.05$ and $p1=.07$, networks from a stochastic block models
with two equal-sized communities and probability link matrices given as

$G\ |\ Y=-1\sim\text{SBM}\begin{pmatrix}.05&.15\\\
.15&.05\end{pmatrix}\quad\text{ and }\quad G\ |\
Y=+1\sim\text{SBM}\begin{pmatrix}.1&.15\\\ .15&.05\end{pmatrix}\enskip,$

and correlated Erdős-Rényi networks (Pedarsani and Grossglauser, 2011) in
which we sample independent Erdős-Rényi networks, $G_{0},G_{1}\sim ER(n,p=.8)$
and then generate $G\ |\ Y=-1$ as

$G_{ij}\ |\ Y=-1\sim\begin{cases}\text{Bernoulli}(\rho=.8)&\text{ if
}G_{0,ij}=1\\\ 0&\text{ if }G_{0,ij}=0\end{cases}\enskip,$

and similarly for $G\ |\ Y=+1$ conditional on $G_{1}$. We also compare
networks from preferential attachment models (Barabási and Albert, 1999) with
different powers $\alpha_{0}=.6$ and $\alpha_{1}=1.4$ and exponential random
graph models that have the same average density and 2-stars but different
transitivity. We refer to these five random graph models as small-world, SBM,
corrER, PA, and ERGM, respectively.

For each network model, we sample $m/2$ networks from each class and use a
75/25 train/test split to assess the classification accuracy. Results are
given as the average accuracy on the test set over 50 replicates in Figure 4,
but the overall conclusions are as follows. For global structure like the
small-world, PA, and ERGM networks, GP-$\lambda$ is the better graph distance,
whereas for local structure like the correlated Erdős-Rényi networks, GP-F is
better. This is unsurprising since the spectral-based distance and Frobenius
distance were expected to capture global and local differences, respectively.
Note that for GP-$\lambda$, the normalized Laplacian consistently outperformed
the Laplacian, so we only report those results and recommend always using the
normalized Laplacian in practice. We also see that the GP-RW is as good or
better than the SVM-RW in all of the graph models. Finally, MASE is an
excellent method in the case when the underlying network structure is low-rank
like the SBM networks. However, we note that we encountered a few instances of
eigendecomposition failure with MASE and those replicates are not included in
the results.

Figure 4: Simulation results for three different random graph models. On the
y-axis, we have classification accuracy using a held out test set. On the
x-axis, we vary the number of networks from $m=40$ to $m=340$, which are split
into two even groups. The column panels correspond to networks with $n$ nodes.

In the first row of Figure 4, we see that GP-$\lambda$ performs the best for
small-world networks across all values of $m$ and $n$. As $n$ increases, most
methods improve with MASE the second most accurate. GP-RW and GC perform
similarly. Note that GP-RW starts to outperform SVM-RW when the ratio $m/n$
exceeds 3-5. This may be due to the fact that we did not scale the grid as the
number of nodes increased, but nevertheless this highlights the flexibility of
our GP model as we use the same hyperprior for all of the settings. Finally,
GP-F is consistently the worst classifier and shows little improvement for
different $m$ and $n$, which is unsurprising since small-world networks are
unlikely to be characterized by local differences.

In the second row, we see that MASE performs the best, though GP-$\lambda$
closes the performance gap as $n$ increases. Similarly, GC is the third best
classifier with performance increasing in both $m$ and $n$. We also see a
clear ranking of GP-RW, then SVM-RW, then GP-F, which again highlights the
superiority of a GP classifier over an SVM despite the similarity of their
objective functions.

The third row shows the results for corrER networks. Note that by design,
these networks are correlated, so this serves as a check of the robustness to
potential violations of independence assumptions. GC and GP-RW rapidly
approach perfect classification accuracy for the different values of $m$ and
$n$. Similar for SVM-RW except for $n=20$ and GP-F except for small $m$.
Interestingly, both MASE and GP-$\lambda$ perform very poorly on this task and
no not improve as $m$ and $n$ increase.

In the fourth and fifth rows, there are no clear trends other than
GP-$\lambda$ outperforming all of the other methods. For the PA networks, all
of the methods improve more with $n$ than with $m$. However, none of the
methods besides GP-$\lambda$ is successful at classifying the ERGM networks,
which appears to be a challenging task.

We also ran our classifiers on the COBRE dataset from Arroyo Relión et al.
(2019). Due to their negative weights, only GP-F is suitable to use these
networks directly. To use GP-$\lambda$, we can represent each network by its
signed Laplacian (Kunegis et al., 2010),

$\bar{L}_{ij}=\begin{cases}-A_{ij}&i\neq j\\\
\sum_{k=1}^{n}|A_{ik}|-A_{ii}&i=j\end{cases}\enskip,$ (8)

and use $\bar{L}$ in (4). Finally, for GP-RW, we need to binarize the network
and we use $0.45$ as a threshold for an edge resulting in $5.3\%$ of the
original edges, which closely matches the best case ($5.4\%$) reported in
Arroyo Relión et al. (2019).

Using the same 10-fold cross-validation as reported by the authors, we
obtained a classification accuracy of $.92$ with a standard deviation of $.10$
for GP-F and also a mean AUC of $.98$ with a standard deviation of $.05$. For
GP-$\lambda$, we obtain an accuracy of $.56$ $(.17)$ and an AUC of $.63$
$(.21)$. Finally, for GP-RW, using the binarized networks, we obtain an
accuracy of $.75$ $(.06)$ and an AUC of $.90$ $(.09)$. These results suggest
that some information may be lost by thresholding the edge weights and also
that the differences between the schizophrenic and the control networks may be
more local than global.

###  5.2 Anomaly Detection

Next, we emulate the scenario when one class is much more prevalent yet we
would like to be able to identify possible samples from another class perhaps
as new samples become available. To do this, we set the classes as unbalanced
(90/10) and train our classifier exclusively on one class, but still use a
75/25 train/test split for evaluating the performance.

From Figure 4, we see that the SBM can be challenging to classify for small
$m$ even when the classes are balanced. We use this regime with $m=100$ and
$n=100$ to illustrate our ability to perform anomaly detection. Note that even
though GC provides classification probabilities, it cannot be trained on a
single class. Likewise for SVM-RW and MASE. Therefore, we only report the
results for GP-F, GP-$\lambda$, and GP-RW, which are given in Figure 5.

Figure 5: Four different OCC scores defined in Section 2.2 applied to our
simulation test set. Scores are derived from the posterior of our GP
classifier trained only on one class (black circle). Dashed lines indicate
possible discrimination thresholds based on the elbow method, i.e. the largest
jump among the sorted values.

We see from the middle row that GP-$\lambda$ has the most discriminatory power
for all of the scores with good separation using the elbow method. GP-F also
has good discrimination power, though the scores are more clumped and the
ordering is reversed. Although GP-RW sorts the scores well, there is little
separation using the elbow method for any of the scores. Overall, the class
probabilities $\tilde{\mu}$ and the coefficient of variation $H$ are the best
scores to use.

###  5.3 Survival Analysis

Finally, we consider two survival analysis problems. In the first case, we
sample $m$ points each from $f_{0}(t)=\mathcal{N}(2,0.8^{2})$ and
$f_{1}(t)=\mathcal{N}(4,1)$, restricted to $\mathbb{R}^{+}$. However, instead
of providing an indicator for which density a point was sampled from, we use
as input a network sampled from $ER(p_{0})$ or $ER(p_{1})$. In this way, we
simultaneously estimate the survival curves and perform classification.
Moreover, we reduce to the noiseless case in the limit as $p_{0}\to 0$ and
$p_{1}\to 1$. In the second case, we consider the more challenging variant
from Fernández, Rivera and Teh (2016), in which
$f_{0}(t)=\mathcal{N}(3,0.8^{2})$ and
$f_{1}(t)=0.4\mathcal{N}(4,1)+0.6\mathcal{N}(2,0.8^{2})$. This is similar to
the first case except that $f_{1}(t)$ is now a mixture of the Gaussians in the
first case, which creates a crossing between the survival functions
corresponding to $f_{0}(t)$ and $f_{1}(t)$.

In both cases, we take $m=50$ samples each from $f_{0}(t)$ and $f_{1}(t)$,
where the network inputs are sampled from an $ER$ with $n=50$ nodes and
$p_{0}=.3$ and $p_{1}=.7$. These classes are very easy to distinguish for both
GP-F and GP-$\lambda$, so we only show the results for GP-F. Also, recall that
our kernel must be stationary for the survival function to be well-defined, so
we cannot use GP-RW.

The results are shown in Figure 6. We see that we successfully estimate
$f_{0}$ and $f_{1}$ in the easy case. On the other hand, the estimates in the
hard case are essentially the mean of $f_{0}$ and $f_{1}$, which are the same
results as De Iorio et al. (2009), who originally introduced this experiment.
Unfortunately, we would need a much larger $m$ to recover the crossing
survival functions, which is not feasible as our algorithm scales poorly in
$m$ even though parallelization helps. Although limiting, we note that our
application only has a total of $m=37$.

Figure 6: Simulation results for the easy case (left) and hard case (right).
The light curves are the posterior means of survival surfaces $S(t,G)$ for
each point, the darker curves are survival curves $S(t)$ of each group and the
black curves are the true survival curves. We sample $m=50$ points from each
density and corresponding ER networks with $n=50$ nodes with $p_{0}=.3$ and
$p_{1}=.7$.

##  6 Application to Microbiome Networks

Finally, we return to our motivating application of the microbiome dataset
from DiGiulio et al. (2015).

###  6.1 Microbiome Network Construction

For every subject in the dataset, there is an operational taxonomic unit (OTU)
table of genomic count data. The rows of an OTU table correspond to the
different samples over the course of a pregnancy and the columns represent
different bacterial taxa. The $(i,j)$-th entry is therefore the count of how
much genomic material of taxa $j$ was found in the sample at time $i$.

Given an OTU table, there are many methods for constructing microbiome
networks. Typically, the nodes of a microbiome network represent a taxa and
the edge weights between nodes represent some association between taxa with
different construction methods reflecting different choices of association.
Herein, the edges of our microbiome networks are correlations, specifically
Sparse Correlations for Compositional data (SparCC) following the well-known
procedure in Friedman and Alm (2012). SparCC approximates the Pearson
correlations between the log-transformed counts with the assumption that the
true correlation network is sparse, though the authors argue their method is
robust to misspecification of the sparsity assumption.

In order to employ the SparCC method, we need to make several preprocessing
decisions on what samples to include. To do so, we follow similar guidelines
as Bogart, Creswell and Gerber (2019), who use the same dataset, among others,
to demonstrate their Bayesian method called MITRE for microbiome time-series
data. In particular, we discard potentially spurious taxa with fewer than 500
total samples and exclude samples where coverage is lower than 10 total
samples. We include only samples from the vaginal swabs and restrict the
samples to between 1 and 30 gestational weeks. While the networks may have
more information if we include samples from all body sites and later time
points, we want to illustrate how they would be used ideally in practice,
which is as a minimally-invasive monitoring procedure throughout pregnancy
that would enable intervention in the case that a preterm delivery may be
likely.

Following this construction, we obtain $m=37$ networks with $n=63$ nodes each.
Two such networks are shown in Figure 2.

###  6.2 Exploratory Analysis

As with studies of non-network data, it is important to first perform an
exploratory analysis. With network data, it is common to try to reduce
networks to a few summary statistics. These summary statistics have been
developed to describe aspects such as the overall network structure, local
connectivity patterns, and important nodes. Although many such statistics have
analogues for weighted networks, the edge weights in our microbiome networks
are correlations and hence can be negative. Unfortunately, methods for
networks with negative weights is an underdeveloped area requiring us to
either use absolute values, which fails to distinguish between negative and
positive associations, or else binarize the weights. Since thresholding is
difficult without having prior information, we use the weight magnitudes for
our exploratory analysis.

In general, the challenge with describing a network is identifying what
statistics are relevant for a given task. In the networks in Figure 2, there
appear to be some visual differences between the preterm and term networks
such as the apparent densities and degree distributions. To assess this, we
compare the averages in the preterm and term groups across several summary
statistics.

We consider four network-level descriptors. We first measure strength, which
is the sum of all of the edge weights in the network. Second, we estimate the
number of communities in each network, which we denote $\hat{K}$. Third, for
distance measures, we use the inverse of the edge weights so that higher
correlations refer to closer nodes. This allows us to measure closeness, a
node-level statistic for how many steps it takes to reach every other node
from a given node, which we aggregate over all nodes. Our final network-level
statistic is diameter, which is the length of the longest path.

We also consider four node-level descriptors. We look at what percentage of
the networks have Gardnerella or Ureaplas, which are both reported in DiGiulio
et al. (2015) to be associated with higher risk of preterm delivery, in their
top three most central nodes. Similarly, we report the percentage for any
Lactobacillus species, which were reported to be inversely correlated with
gestational length, in the top three central nodes. We measure centrality by
eigenvalue, but the results are similar for centrality by strength, closeness,
and hub score. We also report the percentages for Finegoldia, which we have
identified as a central species for term networks that is absent in preterm
networks. This species was not identified previously and deserves further
scientific inquiry. The results are given in Table 1.

Characteristic |  | Term (m = 26) | Preterm (m = 11)
---|---|---|---
Network-level |  |  |
| Strength | 709 (784) | 490 (231)
| $\hat{K}$ | 3.19 (0.69) | 3.45 (0.82)
| Closeness | 0.21 (0.20) | 0.15 (0.06)
| Diameter | 14.0 (6.1) | 15.0 (5.3)
Node-level |  |  |
| Gardnerella | 1/26 (4%) | 3/11 (27%)
| Ureaplas | 2/26 (8%) | 2/11 (18%)
| Lactobacillus | 19/26 (73%) | 2/11 (18%)
| Finegoldia | 5/26 (19%) | 0/11 (0%)
Table 1:  Comparison of average summary statistics for preterm and term
networks. The average (sd) is reported for the network-level statistics and
the proportion of networks with the given species in the top three most
central taxa is reported for the node-level statistics.

Although we see some network-level differences such as strength, none are
significant. This underscores the topological complexity of the microbiome
networks that are clear from Figure 2. We also have some validation of our
network construction by seeing that the keystone species identified in the
literature manifest as differences in node centrality. Moreover, we are also
able to identify a difference in an unreported species. However, without a
priori knowledge of the keystone species, there is an inherent multiple
testing challenge in identifying potentially relevant taxa in new microbiome
networks. Therefore, the entire networks should ideally be used directly as
inputs.

###  6.3 Preterm Delivery Tasks

To apply our models from Section 2, we define three tasks based on the
additional clinical information available in the dataset. In particular, there
is an indicator for a history of preterm delivery, an indicator if the given
pregnancy resulted in a preterm birth, and the length of pregnancy. This
allows us to first define a classification problem using preterm delivery
status. Note that Bogart, Creswell and Gerber (2019) perform the same
classification task, but only group the preterm and very preterm labels in the
dataset as preterm, whereas the original authors included the marginal label
in the preterm class, which they define as before the 37th gestational week.

Next, we note that 11 of the $m=37$ pregnancies in the case-control cohort
resulted in a preterm delivery, but this proportion of preterm deliveries is
higher than the national average of 11%. For this reason, we train our
classifier using a more realistic term/preterm balance and then perform
anomaly detection. We expect this is more likely how our models would be used
in practice i.e given a microbiome network during pregnancy, scores are output
and monitored rather than probabilities. Finally, we perform survival analysis
using days to delivery as our event.

###  6.4 Results

Using leave-one-out cross validation (LOOCV), GP-$\lambda$ using the signed
Laplacian in (8) has a classification accuracy of 78% with an AUC of $0.69$,
whereas GP-F has an accuracy of 61% and an AUC of $0.37$. GP-RW using
binarized networks has an accuracy of 70% with an AUC of $0.53$, but unlike
GP-$\lambda$ and GP-F, it predicts term for every subject. This underscores
the challenge of this classification task. In particular, LOOCV exacerbates
the problem of unbalanced classes, since every time a subject with a preterm
label is withheld, the classes become less balanced. Consequently, metrics
like AUC may not be useful in assessing the classifiers.

Instead, we reformulate the problem as an anomaly detection task by setting
our training set to reflect the population, i.e. roughly 11% of training cases
are preterm. We then run our GP classifier on this unbalanced training set to
obtain OCC scores. The results of this are shown in Figure 7.

Figure 7: OCC scores for the microbiome data. Preterm births are denoted by a
red plus sign and the dashed lines indicate the largest gap in values.

Based on the results, $\tilde{\mu}$ from GP-RW would be most useful. Next, the
low-end scores of $\tilde{\mu}$ and $H$ from GP-F as well as the high-end
scores of $\tilde{\sigma}$ from GP-$\lambda$ If microbiome sampling were to
become part of prenatal care, then these scores could be used to flag high
risk pregnancies rather than using the probabilities of preterm delivery. The
results also suggest that differences in the microbiome networks for preterm
and term pregnancies are a mix of local and global differences, but perhaps
more global due to the slightly higher performance of GP-$\lambda$ over GP-F.
This is consistent with the belief that microbiome networks are complex while
simultaneously supporting the original study findings regarding pronounced
differences in a few taxa.

Finally, we perform survival analysis on time to birth, where we use as input
both the subject’s microbiome network and an indicator for history of preterm
birth. Here, we are able to seamlessly combine data of mixed types using our
graph kernel with the Frobenius distance and a squared-exponential kernel with
Euclidean distance. The result with the spectral-based distance is shown in
Figure 8, which includes a comparison with the Kaplan-Meier (KM) estimates
that only include the indicator as input.

Figure 8: The light curves are the posterior means of survival surfaces
$S(t,G)$ for each subject, the darker curves are survival curves $S(t)$ for
those with and without history of preterm birth, and the dashed curves are
Kaplan-Meier estimates.

One obvious difference between our GP method and the KM estimates is that we
are able to incorporate relevant covariates like the microbiome network and
history of preterm birth. This ability is necessary in practice because
usually there are many other relevant clinical characteristics that need to be
accounted for when personalizing care. Another difference between our GP
method and the KM estimates is that smoothing allows us to extrapolate the
likelihood of delivery prior to any observed time points. This is crucially
important to clinicians because longer gestational periods are directly
related to higher survival rates in pregnancies at higher risk of preterm
delivery.

##  7 Conclusion

###  7.1 Further Applications

In Section 6, we saw how our unified Bayesian framework for classification,
anomaly detection, and survival analysis could be applied to a unique
microbiome dataset. Originally, DiGiulio et al. (2015) demonstrated that
microbiome networks captured the intricate dynamics between microbial taxa in
pregnant women throughout gestation. Our work now represents an important
advance in this research by providing statistical methodology for
incorporating this data as network covariates. Consequently, we hope that this
paper convinces practitioners across various scientific disciplines to collect
network-based data objects and allows them to pose scientific questions using
these networks as covariates. To the best of our knowledge, there are
currently quite few datasets with individual networks and outcomes for all
three problems of classification, anomaly detection, and survival analysis.

We believe that these tools will empower ecologists to design microbiome
studies with many subjects across many time points, enabling further studies
of relationships between the human microbiome and conditions such as pregnancy
or various disease states. Additionally, these tools could be applied to a
variety of non-human microbiomes and soil microbiomes, which are active areas
of research.

Beyond ecology, our work is applicable to other fields that may benefit from
the use of a network per individual or unit. As technological advances are
making this increasingly possible, several areas in particular are poised to
utilize these methods. For example, we have already shown that our models are
useful in neuroscience for classifying diseases using fMRI networks.
Furthermore, these methods can help address clinically relevant time-to-event
questions about the brain such as the development of Alzheimer’s disease.
Other areas of application include the detection of cancer using gene co-
expression networks or single-cell networks, as well as social network
analysis for identifying early adopters or performing time-to-purchase
analyses using ego-nets.

###  7.2 Future Work

There are many exciting extensions to this work. First, it would be impactful
if we had graph distances for networks of different orders or directed
networks. Next, we have seen that kernels seamlessly combine data of mixed
types, such as patient networks and clinical information, but it would be
valuable to have kernels that incorporate exogenous network information.
Similarly, approximations of graph kernels in the spirit of Rahimi and Recht
(2008) would speed up computations and may be necessary as networks become
more present and sample sizes increase. For survival analysis, it may be
possible to adapt the recent variational approach by Kim and Pavlovic (2018)
to network inputs. And, as always, we would benefit from better and faster
MCMC or other Bayesian computational methods.

## References

  * Arroyo et al. (2019) [author] Arroyo, JesúsJ., Athreya, AvantiA., Cape, JoshuaJ., Chen, GuodongG., Priebe, Carey EC. E. and Vogelstein, Joshua TJ. T. (2019). Inference for multiple heterogeneous networks with a common invariant subspace. arXiv preprint arXiv:1906.10026.
  * Arroyo Relión et al. (2019) [author] Arroyo Relión, Jesús DJ. D., Kessler, DanielD., Levina, ElizavetaE., Taylor, Stephan FS. F. et al. (2019). Network classification with applications to brain connectomics. The Annals of Applied Statistics 13 1648–1677.
  * Barabási and Albert (1999) [author] Barabási, Albert-LászlóA.-L. and Albert, RékaR. (1999). Emergence of scaling in random networks. science 286 509–512.
  * Bogart, Creswell and Gerber (2019) [author] Bogart, ElijahE., Creswell, RichardR. and Gerber, Georg KG. K. (2019). MITRE: inferring features from microbiota time-series data linked to host status. Genome biology 20 1–15.
  * Borgwardt et al. (2005) [author] Borgwardt, Karsten MK. M., Ong, Cheng SoonC. S., Schönauer, StefanS., Vishwanathan, SVNS., Smola, Alex JA. J. and Kriegel, Hans-PeterH.-P. (2005). Protein function prediction via graph kernels. Bioinformatics 21 i47–i56.
  * Cornea et al. (2017) [author] Cornea, EmilE., Zhu, HongtuH., Kim, PeterP., Ibrahim, Joseph GJ. G. and Initiative, Alzheimer’s Disease NeuroimagingA. D. N. (2017). Regression models on Riemannian symmetric spaces. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79 463–482.
  * Dai et al. (2018) [author] Dai, XiongtaoX., Müller, Hans-GeorgH.-G. et al. (2018). Principal component analysis for functional data on riemannian manifolds and spheres. The Annals of Statistics 46 3334–3361.
  * De Iorio et al. (2009) [author] De Iorio, MariaM., Johnson, Wesley OW. O., Müller, PeterP. and Rosner, Gary LG. L. (2009). Bayesian nonparametric nonproportional hazards survival modeling. Biometrics 65 762–771.
  * DiGiulio et al. (2015) [author] DiGiulio, Daniel BD. B., Callahan, Benjamin JB. J., McMurdie, Paul JP. J., Costello, Elizabeth KE. K., Lyell, Deirdre JD. J., Robaczewska, AnnaA., Sun, Christine LC. L., Goltsman, Daniela SAD. S., Wong, Ronald JR. J., Shaw, GaryG. et al. (2015). Temporal and spatial variation of the human microbiota during pregnancy. Proceedings of the National Academy of Sciences 112 11060–11065.
  * Donnat et al. (2018) [author] Donnat, ClaireC., Holmes, SusanS. et al. (2018). Tracking network dynamics: A survey using graph distances. The Annals of Applied Statistics 12 971–1012.
  * Durante, Dunson and Vogelstein (2017) [author] Durante, DanieleD., Dunson, David BD. B. and Vogelstein, Joshua TJ. T. (2017). Nonparametric Bayes modeling of populations of networks. Journal of the American Statistical Association 112 1516–1530.
  * Fernández, Rivera and Teh (2016) [author] Fernández, TamaraT., Rivera, NicolásN. and Teh, Yee WhyeY. W. (2016). Gaussian processes for survival analysis. In Advances in Neural Information Processing Systems 5021–5029.
  * Fernández and Teh (2016) [author] Fernández, TamaraT. and Teh, Yee WhyeY. W. (2016). Posterior Consistency for a Non-parametric Survival Model under a Gaussian Process Prior. arXiv preprint arXiv:1611.02335.
  * Friedman and Alm (2012) [author] Friedman, JonathanJ. and Alm, Eric JE. J. (2012). Inferring correlation networks from genomic survey data. PLoS computational biology 8\.
  * Gärtner, Driessens and Ramon (2003) [author] Gärtner, ThomasT., Driessens, KurtK. and Ramon, JanJ. (2003). Graph kernels and gaussian processes for relational reinforcement learning. In International Conference on Inductive Logic Programming 146–163. Springer.
  * Ghosal and Roy (2006) [author] Ghosal, SubhashisS. and Roy, AnindyaA. (2006). Posterior consistency of Gaussian process prior for nonparametric binary regression. The Annals of Statistics 2413–2429.
  * Ginestet et al. (2017) [author] Ginestet, Cedric EC. E., Li, JunJ., Balachandran, PrakashP., Rosenberg, StevenS., Kolaczyk, Eric DE. D. et al. (2017). Hypothesis testing for network data in functional neuroimaging. The Annals of Applied Statistics 11 725–750.
  * Gollini and Murphy (2016) [author] Gollini, IsabellaI. and Murphy, Thomas BrendanT. B. (2016). Joint modeling of multiple network views. Journal of Computational and Graphical Statistics 25 246–265.
  * Hsu et al. (2003) [author] Hsu, Chih-WeiC.-W., Chang, Chih-ChungC.-C., Lin, Chih-JenC.-J. et al. (2003). A practical guide to support vector classification.
  * Jain (2016) [author] Jain, Brijnesh JB. J. (2016). On the geometry of graph spaces. Discrete Applied Mathematics 214 126–144.
  * Jayasumana et al. (2013) [author] Jayasumana, SadeepS., Hartley, RichardR., Salzmann, MathieuM., Li, HongdongH. and Harandi, MehrtashM. (2013). Kernel methods on the Riemannian manifold of symmetric positive definite matrices. In Computer Vision and Pattern Recognition (CVPR), 2013 IEEE Conference on 73–80. IEEE.
  * Kashima and Inokuchi (2002) [author] Kashima, HisashiH. and Inokuchi, AkihiroA. (2002). Kernels for graph classification. In ICDM workshop on active mining 2002\.
  * Kemmler et al. (2013) [author] Kemmler, MichaelM., Rodner, ErikE., Wacker, Esther-SabrinaE.-S. and Denzler, JoachimJ. (2013). One-class classification with Gaussian processes. Pattern Recognition 46 3507–3518.
  * Khan and Madden (2009) [author] Khan, Shehroz SS. S. and Madden, Michael GM. G. (2009). A survey of recent trends in one class classification. In Irish conference on artificial intelligence and cognitive science 188–197. Springer.
  * Kim and Pavlovic (2018) [author] Kim, MinyoungM. and Pavlovic, VladimirV. (2018). Variational Inference for Gaussian Process Models for Survival Analysis. In UAI 435–445.
  * Kolaczyk and Csárdi (2014) [author] Kolaczyk, Eric DE. D. and Csárdi, GáborG. (2014). Statistical analysis of network data with R 65\. Springer.
  * Kolaczyk et al. (2020) [author] Kolaczyk, Eric DE. D., Lin, LizhenL., Rosenberg, StevenS., Walters, JacksonJ. and Xu, JieJ. (2020). Averages of unlabeled networks: Geometric characterization and asymptotic behavior. The Annals of Statistics 48 514–538.
  * Kondor and Lafferty (2002) [author] Kondor, Risi ImreR. I. and Lafferty, JohnJ. (2002). Diffusion kernels on graphs and other discrete structures. In Proceedings of the 19th international conference on machine learning 2002 315–322.
  * Kriege, Johansson and Morris (2020) [author] Kriege, Nils MN. M., Johansson, Fredrik DF. D. and Morris, ChristopherC. (2020). A survey on graph kernels. Applied Network Science 5 1–42.
  * Kunegis et al. (2010) [author] Kunegis, JérômeJ., Schmidt, StephanS., Lommatzsch, AndreasA., Lerner, JürgenJ., De Luca, Ernesto WE. W. and Albayrak, SahinS. (2010). Spectral analysis of signed graphs for clustering, prediction and visualization. In Proceedings of the 2010 SIAM International Conference on Data Mining 559–570. SIAM.
  * Layeghifard, Hwang and Guttman (2017) [author] Layeghifard, MehdiM., Hwang, David MD. M. and Guttman, David SD. S. (2017). Disentangling interactions in the microbiome: a network perspective. Trends in microbiology 25 217–228.
  * McNeish (2016) [author] McNeish, DanielD. (2016). On using Bayesian methods to address small sample problems. Structural Equation Modeling: A Multidisciplinary Journal 23 750–773.
  * Mukherjee, Sarkar and Lin (2017) [author] Mukherjee, Soumendu SundarS. S., Sarkar, PurnamritaP. and Lin, LizhenL. (2017). On clustering network-valued data. In Advances in Neural Information Processing Systems (I.I. Guyon, U. V.U. V. Luxburg, S.S. Bengio, H.H. Wallach, R.R. Fergus, S.S. Vishwanathan and R.R. Garnett, eds.) 30 7071–7081.
  * Murray and Adams (2010) [author] Murray, IainI. and Adams, Ryan PR. P. (2010). Slice sampling covariance hyperparameters of latent Gaussian models. In Advances in neural information processing systems 1732–1740.
  * Murray, Adams and MacKay (2010) [author] Murray, IainI., Adams, Ryan PR. P. and MacKay, David JCD. J. (2010). Elliptical slice sampling.
  * Mygdalis et al. (2016) [author] Mygdalis, VasileiosV., Iosifidis, AlexandrosA., Tefas, AnastasiosA. and Pitas, IoannisI. (2016). Graph embedded one-class classifiers for media data classification. Pattern Recognition 60 585–595.
  * Nikolentzos, Siglidis and Vazirgiannis (2019) [author] Nikolentzos, GiannisG., Siglidis, GiannisG. and Vazirgiannis, MichalisM. (2019). Graph Kernels: A Survey. arXiv preprint arXiv:1904.12218.
  * Pedarsani and Grossglauser (2011) [author] Pedarsani, PedramP. and Grossglauser, MatthiasM. (2011). On the privacy of anonymized networks. In Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining 1235–1243.
  * Rahimi and Recht (2008) [author] Rahimi, AliA. and Recht, BenjaminB. (2008). Random features for large-scale kernel machines. In Advances in neural information processing systems 1177–1184.
  * Ralaivola et al. (2005) [author] Ralaivola, LivaL., Swamidass, Sanjay JS. J., Saigo, HirotoH. and Baldi, PierreP. (2005). Graph kernels for chemical informatics. Neural networks 18 1093–1110.
  * Rasmussen and Williams (2005) [author] Rasmussen, Carl EdwardC. E. and Williams, Christopher K. I.C. K. I. (2005). Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press.
  * Rudd (2018) [author] Rudd, Jessica MJ. M. (2018). Application of support vector machine modeling and graph theory metrics for disease classification. Model Assisted Statistics and Applications 13 341–349.
  * Salter-Townshend and McCormick (2017) [author] Salter-Townshend, MichaelM. and McCormick, Tyler HT. H. (2017). Latent space models for multiview network data. The annals of applied statistics 11 1217\.
  * Tang et al. (2018) [author] Tang, RunzeR., Ketcha, MichaelM., Badea, AlexandraA., Calabrese, Evan DE. D., Margulies, Daniel SD. S., Vogelstein, Joshua TJ. T., Priebe, Carey EC. E. and Sussman, Daniel LD. L. (2018). Connectome smoothing via low-rank approximations. IEEE transactions on medical imaging 38 1446–1456.
  * Tokdar and Ghosh (2007) [author] Tokdar, Surya TS. T. and Ghosh, Jayanta KJ. K. (2007). Posterior consistency of logistic Gaussian process priors in density estimation. Journal of statistical planning and inference 137 34–42.
  * Vishwanathan et al. (2010) [author] Vishwanathan, S Vichy NS. V. N., Schraudolph, Nicol NN. N., Kondor, RisiR. and Borgwardt, Karsten MK. M. (2010). Graph kernels. Journal of Machine Learning Research 11 1201–1242.
  * Watts and Strogatz (1998) [author] Watts, Duncan JD. J. and Strogatz, Steven HS. H. (1998). Collective dynamics of ‘small-world’networks. nature 393 440–442.
  * Zhang et al. (2013) [author] Zhang, WeiW., Ota, TakayoT., Shridhar, VijiV., Chien, JeremyJ., Wu, BaolinB. and Kuang, RuiR. (2013). Network-based survival analysis reveals subnetwork signatures for predicting outcomes of ovarian cancer treatment. PLoS computational biology 9 e1002975.

##

###### Proof of Theorem 4.2.

Let $\varepsilon>0$. We want to show that

$\Pi\Big{(}p:\int|p(x)-p_{0}(x)|dQ(x)>\varepsilon\Big{|}Y_{1},\ldots,Y_{m},X_{1},\ldots,X_{m}\Big{)}\to
0\enskip,$

which, by Ghosal and Roy (2006), happens just in case the following four
conditions hold:

  * (P)

For every $x\in\mathcal{X}$, the covariance function $\sigma_{0}(x,\cdot)$ has
continuous partial derivatives up to order $2\alpha+2$, the mean function
$\mu(x)\in\bar{\mathcal{A}}$, the RKHS of $\sigma_{0}(\cdot,\cdot)$, and the
prior $\Pi_{\lambda}$ is fully supported on $(0,\infty)$ ;

  * (C)

The covariate space $\mathcal{X}$ is a bounded subset of $\mathbb{R}^{d}$ ;

  * (T)

The transformed true response function $\eta_{0}\in\bar{\mathcal{A}}$ ;

  * (G)

For every $b_{1}>0$ and $b_{2}>0$, there exist sequences $M_{n},\tau_{n}$ and
$\lambda_{n}$ such that

$M_{m}^{2}\tau_{m}\lambda_{m}^{-2}\geq b_{1}m\qquad\text{ and }\qquad
M_{m}^{d/\alpha}\leq b_{2}m\enskip.$

It is straightforward to verify conditions $(P)$, $(C)$, and $(G)$:

  * (P)

The squared-exponential kernel is infinitely divisible, so we can take any
$\alpha\in\mathbb{N}$ in $(P)$. Furthermore, we are taking
$\mu(x)=0\in\bar{\mathcal{A}}$ since all RKHS contain the identity element.
Finally, we put an inverse-gamma prior on the bandwidth $\ell$, which is fully
supported on the positive reals;

  * (C)

The boundedness of our covariate space follows from the fact that we can embed
the space of weighted networks in $\mathbb{R}^{n^{2}}$, where $n$ is the
number of nodes;

  * (G)

We can take $M_{m}=\ell_{m}$ and $\tau_{m}=b_{1}m$. Then the result holds,
since we can take any $\alpha\in\mathbb{N}$ .

Note that the embedding of the space of weighted networks in
$\mathbb{R}^{n^{2}}$ is given by Ginestet et al. (2017).

Our only remaining obstacle to showing posterior consistency is $(T)$. Tokdar
and Ghosh (2007) provide an equivalency for $(T)$. The authors show that

$\eta_{0}\in\bar{\mathcal{A}}\iff\forall\varepsilon>0,\
\mathbb{P}(\|\eta(x)-\eta_{0}\|_{\infty}<\varepsilon)>0$

if the following three conditions hold:

  * (A1)

$\exists M,M^{\prime}>0$ such that $M^{\prime}\leq\sigma_{0}(t,t)\leq M,\
\forall t\in(\mathbb{R}^{+})^{d}$ ;

  * (A2)

$\exists C,q>0$ such that
$[\sigma_{0}(t,t)+\sigma_{0}(s,s)-2\sigma_{0}(t,s)]^{1/2}\leq C\|s-t\|^{q},\
\forall s,t\in(\mathbb{R}^{+})^{d}$ ;

  * (A3)

For any $n\geq 1$ and any
$t_{1},\ldots,t_{n}\in(\mathbb{R}^{+})^{d},\Sigma=((\sigma(t_{i},t_{j})))$ is
nonsingular.

First, we verify these conditions.

  * (A1)

Again, since we are using the squared-exponential kernel, for any distance
$d(\cdot,\cdot)$, we have $d(t,t)=0$ so that
$\exp(-d^{2}(t,t)/2\sigma^{2})=1$. Therefore, take $M^{\prime}=1=M$ ;

  * (A2)

Since $\sigma_{0}(t,s)\leq 1$, we can take $q=1$ to obtain

$[1+1-2\sigma_{0}(t,s)]^{1/2}\leq\sqrt{2}\enskip.$

Therefore, take $C=\sqrt{2}/\max_{s,t\in(\mathbb{R}^{+})^{d}}|s-t|$ ;

  * (A3)

We already proved that a squared-exponential kernel induced by the Frobenius
distance in (3) or the spectral-based distance in (4) is positive-definite.

So, we are done if we can show

$\forall\varepsilon>0,\
\mathbb{P}(\|\eta(x)-\eta_{0}\|_{\infty}<\varepsilon)>0\enskip.$

For $\varepsilon>0$, we have

$\displaystyle\mathbb{P}(\|\eta(x)-\eta_{0}\|_{\infty}<\varepsilon)$
$\displaystyle\geq\mathbb{P}(\|\eta(x)\|_{\infty}+\|\eta_{0}\|_{\infty}<\varepsilon)=\mathbb{P}(\|\eta(x)\|_{\infty}<\varepsilon-|\eta_{0}|)$
$\displaystyle=\mathbb{P}(\sup_{x}|\eta(x)|<\varepsilon-|\eta_{0}|)=1-\mathbb{P}(\sup_{x}|\eta(x)|>\varepsilon-|\eta_{0}|)$
$\displaystyle\geq 1-\exp(-(\varepsilon-|\eta_{0}|)^{2}/2)>0\enskip,$

where we used the Borell, or Borell-TIS, inequality in the fifth step, which
says that for a mean-zero GP, $X$, with $\sigma^{2}(X))=\sup\text{Var}(X)$, we
have

$\mathbb{P}(\|X\|_{\infty}\geq x)\leq 2\exp(-x^{2}/2\sigma^{2}(X))\enskip.$

Hence, $\eta_{0}\in\bar{\mathcal{A}}$, verifying condition $(T)$ and
concluding our proof. ∎

###### Proof of Theorem 4.3.

Fernández and Teh (2016) show that for a randomized design, the model for
survival analysis in Section 2.3 achieves posterior consistency if the
following four conditions hold:

  * (A1)

$\big{(}k(0)-k(2^{-n})\big{)}^{-1}\geq n^{6}$ ;

  * (A2)

$\nu$ assigns positive probability to every neighborhood of the true
parameter;

  * (A3)

$\mathbb{E}_{\theta_{0}}[T]<M<\infty$ ;

  * (A4)

The true parameters $\eta_{j,0}$ take the form $\hat{\eta}_{i,0}/h_{d}$, where
$\hat{\eta}_{j,0}$ is in the support of $\hat{\eta}_{i}$ under the uniform
norm. $\hat{\eta}_{j}=h_{d}\eta_{j}$ is the non-stationary GP as a function of
the original GP with

$h_{d}(t)=\begin{cases}\frac{d+1}{1+\log(1-e^{-1})}&t<1\\\
\frac{d+1}{t+\log(1-e^{-t})}&t\geq 1\enskip.\end{cases}$

The authors state that these assumptions are “quite reasonable and natural for
the type of data we were dealing with.” In particular, we assume (A3) and (A4)
hold. In words, (A3) just says that survival times have finite expectation,
which is true for survival times like death. For others, like diagnosis, it is
ill-posed to consider whether the likelihood goes to one as time increases. In
general, this is a very reasonable assumption. On the other hand, (A4) is
unverifiable, so we take it on faith.

For (A1), we note that for squared-exponential kernels, which are stationary,
we have

$k(s)=\sigma^{2}\cdot\exp\Big{(}-\frac{s^{2}}{2\ell^{2}}\Big{)}=\sigma^{2}\cdot\exp\Big{(}-\tilde{\ell}\cdot
s^{2}\Big{)}\enskip,$

for some distance $s$. Hence, there is no ambiguity about what stationarity
means for graph kernels, as the stationarity is with respect to, in our case,
both the Frobenius distance and the spectral-based distance. Therefore, we
have

$\big{(}k(0)-k(2^{-n})\big{)}^{-1}\geq
n^{6}\iff\frac{\big{(}1-\exp(-\frac{\ell}{4^{n}})\big{)}^{-1}}{\sigma^{2}\cdot
n^{6}}\geq 1\enskip.$

This is true asymptotically, as the limit of the left side of the inequality
goes to $\infty$ in $n$.

For (A2) We take a constant baseline hazard function $\omega=2\cdot\Omega$.
This corresponds to a hazard function of an exponential random variable with
mean $1/\Omega$, i.e. $\omega_{0}\sim\nu=\text{Exp}(\Omega)$, which is fully
supported on the positive reals. ∎