#  ssROC: Semi-Supervised ROC Analysis for
Reliable and Streamlined Evaluation of Phenotyping Algorithms

Jianhui Gao† Department of Statistical Sciences, University of Toronto,
Toronto, ON, Canada  Clara-Lea Bonzel† Department of Biomedical Informatics,
Harvard Medical School, Boston, MA, USA  Chuan Hong  Department of
Biostatistics and Bioinformatics, Duke University, Durham, NC,USA  Paul
Varghese  Verily Life Sciences, Cambridge, MA, USA  Karim Zakir  Department of
Statistical Sciences, University of Toronto, Toronto, ON, Canada  Jessica
Gronsbell  Department of Statistical Sciences, University of Toronto, Toronto,
ON, Canada  Department of Family and Community Medicine, University of
Toronto, Toronto, ON, Canada  Department of Computer Science, University of
Toronto, Toronto, ON, Canada

Correspondence to:
Jessica Gronsbell
Postal address: 700 University Ave, Toronto, ON, Canada, M5G 1Z5
Email<EMAIL_ADDRESS>
Telephone number: 416-978-3452

Keywords: Electronic Health Records; Phenotyping; Semi-supervised; ROC
Analysis
Word count: 3990/4000

## ABSTRACT

Objective: High-throughput phenotyping will accelerate the use of electronic
health records (EHRs) for translational research. A critical roadblock is the
extensive medical supervision required for phenotyping algorithm (PA)
estimation and evaluation. To address this challenge, numerous weakly-
supervised learning methods have been proposed to estimate PAs. However, there
is a paucity of methods for reliably evaluating the predictive performance of
PAs when a very small proportion of the data is labeled. To fill this gap, we
introduce a semi-supervised approach (ssROC) for estimation of the receiver
operating characteristic (ROC) parameters of PAs (e.g., sensitivity,
specificity).

Materials and methods: ssROC uses a small labeled dataset to nonparametrically
impute missing labels. The imputations are then used for ROC parameter
estimation to yield more precise estimates of PA performance relative to
classical supervised ROC analysis (supROC) using only labeled data. We
evaluated ssROC through in-depth simulation studies and an extensive
evaluation of eight PAs from Mass General Brigham.

Results: In both simulated and real data, ssROC produced ROC parameter
estimates with significantly lower variance than supROC for a given amount of
labeled data. For the eight PAs, our results illustrate that ssROC achieves
similar precision to supROC, but with approximately 60% of the amount of
labeled data on average.

Discussion: ssROC enables precise evaluation of PA performance to increase
trust in observational health research without demanding large volumes of
labeled data. ssROC is also easily implementable in open-source R software.

Conclusion: When used in conjunction with weakly-supervised PAs, ssROC
facilitates the reliable and streamlined phenotyping necessary for EHR-based
research.

## BACKGROUND AND SIGNIFICANCE

Electronic Health Records (EHRs) are a vital source of data for clinical and
translational research [1]. Vast amounts of patient data are being tapped for
real-time studies of infectious diseases, development of clinical decision
support tools, and genetic studies at unprecedented scale [2, 3, 4, 5, 6, 7,
8, 9, 10]. This myriad of opportunities rests on the ability to accurately and
rapidly extract a wide variety of patient phenotypes (e.g., diseases) to both
identify and characterize populations of interest. However, precise and
readily available phenotype information is rarely available in patient records
and presents a major barrier to EHR-based research [11, 12].

In practice, phenotypes are extracted from patient records with either rule-
based or machine learning (ML)-based phenotyping algorithms (PAs) derived from
codified and natural language processing (NLP)-derived features [13, 14].
While PAs can characterize clinical conditions with high accuracy, they
traditionally require a substantial amount of medical supervision that limits
the automated power of EHR-based studies [15]. Several research networks have
spent considerable effort developing PAs over the last decade, including i2b2
(Informatics for Integrating Biology & the Bedside), the eMERGE (Electronic
Medical Records and Genomics) Network, and the OHDSI (Observational Health
Data Sciences and Informatics) program that released the APHRODITE (Automated
PHenotype Routine for Observational Definition, Identification, Training, and
Evaluation) framework [16].

Typically, PA development consists of two key steps: (i) algorithm estimation
and (ii) algorithm evaluation. Algorithm estimation determines the appropriate
aggregation of features extracted from patient records to determine phenotype
status. For a rule-based approach, medical and informatics experts assemble a
comprehensive set of features and corresponding logic to assign patient
phenotypes [12, 13]. As this manual assembly can be highly laborious,
significant effort has been made to automate the algorithm estimation step
with ML. Numerous studies have demonstrated success with PAs derived from
standard supervised learning methods such as penalized logistic regression,
random forest, and deep neural networks [17, 18, 19, 20, 21, 22, 23, 24]. The
scalability of a supervised approach, however, is limited by the substantial
number of gold-standard labels required for model training. Gold-standard
labels, which require time-consuming manual medical chart review, are
infeasible to obtain for a large volume of records [25, 26].

In response, semi-supervised and weakly-supervised methods for PA estimation
that substantially decrease or eliminate the need for gold-standard labeled
data have been proposed. Among semi-supervised methods, self-learning and
surrogate-assisted semi-supervised learning are common [27, 28, 29]. For
example, [15] introduced PheCAP, a common pipeline for semi-supervised
learning that utilizes silver-standard labels for feature selection prior to
supervised model training to decrease the labeling demand. Unlike gold-
standard labels, silver-standard labels can be automatically extracted from
patient records (e.g., ICD codes or free-text mentions of the disease) and
serve as proxies for the gold-standard label [30, 31]. PheCAP was based on the
pioneering work of [32] and [33], which introduced weakly-supervised PAs
trained entirely on silver-standard labels. These methods completely eliminate
the need for medical chart review for algorithm estimation and are the basis
of the APHRODITE framework. Moreover, this work prompted numerous developments
in weakly-supervised PAs, including methods based on non-negative
matrix/tensor factorization, parametric mixture modeling, and neural network-
based embedding learning, which are quickly becoming the new standard in the
PA literature [14].

In contrast to the success in automating PA estimation, there has been little
focus on the algorithm evaluation step. Algorithm evaluation assesses the
predictive performance of a PA, typically through the estimation of the
receiver operating characteristic (ROC) parameters such as sensitivity and
specificity. At a high-level, the ROC parameters measure how well a PA
discriminates between phenotype cases and controls relative to the gold-
standard obtained from chart review. As phenotypes are the foundation of EHR-
based studies, it is critical to reliably evaluate the ROC parameters to
provide researchers with a sense of trust in using a PA [34, 35, 36]. However,
complete PA evaluation is performed far too infrequently due to the burden of
chart review [25, 37, 14].

To address this challenge, [38] proposed the first semi-supervised method for
ROC parameter estimation. This method assumes that the predictive model is
derived from a penalized logistic regression model and was only validated on 2
PAs with relatively large labeled data sets (455 and 500 labels). [25] later
introduced PheValuator, and its recent successor PheValuator 2.0, to
efficiently evaluate rule-based algorithms using “probabilistic gold-standard”
labels generated from diagnostic predictive models rather than chart review
[37]. Although the authors provided a comprehensive evaluation for numerous
rule-based PAs, PheValuator can lead to biased ROC analysis, and hence a
distorted understanding of the performance of a PA, when the diagnostic
predictive model is not well-calibrated [39, 40, 41]. PheValuator is also only
applicable for the evaluation of rule-based PAs.

To fill this gap in the PA literature, we introduce the semi-supervised
approach of [38] to precisely estimate the ROC parameters of PAs, which we
call “ssROC”. The key difference between ssROC and classical ROC analysis
(supROC) using only labeled data is that ssROC imputes missing gold-standard
labels in order to leverage large volumes of unlabeled data (i.e., records
without gold-standard labels). By doing so, ssROC yields less variable
estimates than supROC to enable reliable PA evaluation with fewer gold-
standard labels and hence significantly less chart review. Moreover, ssROC
imputes the missing labels with a nonparametric recalibration of the
predictions from the PA to ensure that the resulting estimates of the ROC
parameters are unbiased when the PA is not well-calibrated.

## OBJECTIVE

The primary objectives of this work are to:

  1. 1.

Extend the proposal of [38] to a wider class of weakly-supervised PAs that are
commonly used in the PA literature, including a theoretical analysis and
development of a statistical inference procedure demonstrating good finite-
sample properties.

  2. 2.

Provide an in-depth real data analysis of PAs for eight phenotypes from Mass
General Brigham (MGB) together with extensive simulation studies to illustrate
the practical utility of ssROC.

  3. 3.

Release an implementation of ssROC in open-source R software to encourage the
use of our method by the informatics community.

Through our real and simulated data analyses, we observe substantial gains in
estimation precision from ssROC relative to supROC. In the analysis of the
eight PAs from MGB, ssROC achieves similar precision to supROC, but with
approximately 60% of the amount of labeled data on average. Our results
suggest that, when used together with weakly supervised PAs, ssROC can
facilitate the reliable and streamlined PA development that is necessary for
EHR-based research.

## MATERIALS AND METHODS

### Overview of ssROC

We focus on evaluating a classification rule derived from a PA with ROC
analysis. ROC analysis assesses the agreement between the gold-standard label
for a binary phenotype (e.g., disease case/control), $Y$, and a PA score, $S$,
indicating a patient’s likelihood of having the underlying phenotype (e.g.,
the predicted probability of being a disease case). $Y$ is typically obtained
from chart review and $S$ can be derived from various phenotyping methods. In
this article, we primarily consider PAs derived from weakly-supervised
learning methods due to their ability to automate PA estimation and increasing
popularity in the informatics literature [14].

In classical supervised ROC analysis, it is assumed that the data is fully
labeled and one has information on both $Y$ and $S$ for all subjects. However,
in the phenotyping setting, $Y$ is typically only available for a very small
subset of patient records due to the laborious nature of medical chart review.
This gives rise to the semi-supervised setting in which a small labeled
dataset is accompanied by a much larger unlabeled dataset. To fully leverage
all of the available data and facilitate more reliable (i.e., lower variance)
evaluation of PAs, ssROC imputes the missing $Y$ with a nonparametric
recalibration of $S$, denoted as ${\widehat{m}}(S)$, to make use of the
unlabeled data. A schematic representation of ssROC is provided in Figure 1.

Figure 1: Overview of PA estimation and evaluation. The phenotyping algorithm
(PA) is first estimated to obtain the scores ($S$). Patient charts from the
electronic health record (EHR) warehouse are also reviewed to obtain the gold-
standard label ($Y$) for PA evaluation. In classical supervised ROC analysis
(supROC), only the labeled data from chart review is used to evaluate the PA’s
performance. Semi-supervised ROC analysis (ssROC), uses the labeled data to
impute the missing $Y$ as ${\widehat{m}}(S)$ so that the unlabeled data can be
utilized for estimation to yield more precise estimates of the ROC parameters.

### Data structure & notation

More concretely, the available data in the semi-supervised setting consists of
a small labeled dataset with independently and identically distributed
observations,

$\mathscr{L}=\left\\{(Y_{i},S_{i})\mid i=1,\ldots,n\right\\}$

and an independent unlabeled dataset

$\mathscr{U}=\left\\{S_{i}\mid i=n+1,\ldots,n+N\right\\}.$

Additionally, (i) $\mathscr{U}$ is assumed to be much larger than
$\mathscr{L}$ so that $n\ll N$ and (ii) the observations in $\mathscr{L}$ are
randomly selected from the underlying pool of data. Throughout our discussion,
we suppose that a higher value of $S$ is more indicative of the phenotype. An
observation is deemed to have the phenotype if $S>c$, where $c$ is the
threshold used for phenotype classification.

### ROC analysis

More formally, ROC analysis evaluates a PA based on the true positive rate
(TPR), the false positive rate (FPR), the positive predictive value (PPV), and
the negative predictive value (NPV). In the diagnostic testing literature, the
TPR is referred to as sensitivity, while the FPR is 1 minus the specificity
[42]. For a given classification threshold, one may evaluate the ROC
parameters by enumerating the correct and incorrect classifications. This
information can be summarized in a confusion matrix, as in Figure 2.

Figure 2: Confusion matrix. The algorithm score from the PA is used to
determine phenotype case/control status based on the classification threshold.
The ROC parameters can then be evaluated by enumerating the number of correct
and incorrect classifications relative to the gold-standard label.

In practice, it is the task of the researcher to estimate an appropriate
threshold for classification. This is commonly done by summarizing the trade-
off between the TPR and FPR, defined as

$\displaystyle\textrm{TPR}(c)=P(S>c\ \mid\ Y=1)\text{ and
}\textrm{FPR}(c)=P(S>c\ \mid\ Y=0).$

The ROC curve, $\textrm{ROC}(u)=\textrm{TPR}[\textrm{FPR}^{-1}(u)]$,
summarizes the TPR and FPR across all possible choices of the threshold. In
the context of PAs, $c$ is often chosen to achieve a low FPR, particularly
when the prevalence of the phenotype is low [22]. An overall summary measure
of the discriminative power of $S$ in classifying $Y$ is captured by the area
under the ROC curve (AUC),

$\mbox{AUC}=\int_{0}^{1}\textrm{ROC}(u)du.$

The AUC is equivalent to the probability that a phenotype case has a higher
value of $S$ than a phenotype control. For a given threshold, the predictive
performance of the classification rule derived from the PA is assessed with
the PPV and NPV, defined as

$\displaystyle\textrm{PPV}(c)=P(Y=1\ \mid\ S>c),$
$\displaystyle\quad\mbox{and}\quad\textrm{NPV}(c)=P(Y=0\ \mid\ S<c).$

### supROC: classical supervised ROC analysis

With only labeled data, one may obtain supervised estimators of the ROC
parameters (supROC) based on their empirical counterparts. For example, the
TPR and FPR can be estimated as

$\displaystyle\widehat{\textrm{TPR}}_{sup}(c)=\frac{\sum_{i=1}^{n}Y_{i}I(S_{i}>c)}{\sum_{i=1}^{n}Y_{i}}\text{
and
}\widehat{\textrm{FPR}}_{sup}(c)=\frac{\sum_{i=1}^{n}(1-Y_{i})I(S_{i}>c)}{\sum_{i=1}^{n}(1-Y_{i})}$

for a given $c$. The PPV and NPV can be estimated analogously. A point along
the ROC curve is given by
$\widehat{\textrm{ROC}}_{sup}(u_{0})=\widehat{\textrm{TPR}}_{sup}(\hat{c}_{u_{0}})=\widehat{\textrm{TPR}}_{sup}\left[\widehat{\textrm{FPR}}_{sup}^{-1}(u_{0})\right]$
while the AUC may be obtained from standard numerical integration of the
estimated ROC curve [43]. Variance estimates are obtained with a resampling-
based procedure such as bootstrap or perturbation resampling [44, 45].

### ssROC: semi-supervised ROC Analysis

Unlike its supervised counterpart that relies on only $\mathscr{L}$, ssROC
utilizes both $\mathscr{L}$ and $\mathscr{U}$ to yield more precise estimates
of the ROC parameters. To do so, ssROC consists of two steps of estimation
(Figure 1). In the first step, the missing labels in $\mathscr{U}$ are imputed
with a nonparametric recalibration of the PA scores via local constant
regression with $\mathscr{L}$. In contrast to using the PA scores directly for
imputation, this step ensures that the ROC estimates are robust to potential
bias from poor calibration of the PA [40]. Intuitively, this choice of
imputation provides ssROC with the best guess of $Y$ based on $S$ for ROC
analysis [40]. In the second step, ssROC uses the imputations and PA scores in
$\mathscr{U}$ to estimate the ROC parameters. This step produces estimates
that have lower variance than supROC as $\mathscr{U}$ is assumed to be much
larger than $\mathscr{L}$. Practically, the purpose of ssROC is to enable a
more reliable understanding of a PA’s performance by yielding less variable
estimates relative to supROC without requiring more labeled data. Using the
TPR for illustration, the ssROC method is detailed in Figure 3.

Overview of ssROC Step 1. Estimate $m(s)=P(Y=1\mid S=s)$ with $\mathscr{L}$
via local constant regression as
$\widehat{m}(s)=\frac{\sum_{i=1}^{n}K_{h}(S_{i}-s)Y_{i}}{\sum_{i=1}^{n}K_{h}(S_{i}-s)}$
where $K_{h}(u)=h^{-1}K(u/h)$ is a given smooth, symmetric kernel function,
$h$ is a bandwidth controlling the smoothing, and $nh^{2}\to\infty$ and
$nh^{4}\to 0$ as $n\to\infty$. Step 2. Impute the missing labels in
$\mathscr{U}$ with ${\widehat{m}}(\cdot)$ and evaluate the TPR analogously to
the supervised setting with the imputations in place of the gold-standard
labels as
$\widetilde{\textrm{TPR}}_{ssROC}(c)=\frac{\sum_{i=n+1}^{n+N}\widehat{m}(S_{i})I(S_{i}>c)}{\sum_{i=n+1}^{n+N}\widehat{m}(S_{i})}.$
Figure 3: Two steps of ssROC using the true positive rate (TPR) as an example.

The conditions on the bandwidth in Step 1 ensure that $\widehat{m}(\cdot)$ is
undersmoothed to prevent bias in the estimated ROC parameter [46].
Additionally, the imputations obtained in Step 1 can be used in an analogous
manner to compute the remaining ROC parameters. Similar to supROC analysis, we
propose a perturbation resampling procedure for variance estimation and detail
two commonly used confidence intervals (CIs) based on the resampling procedure
in Supplementary Section 1. In Supplementary Section 2, we also provide a
brief theoretical justification for the improved efficiency of ssROC relative
to supROC for a wide variety of weakly-supervised PAs. Our analysis verifies
that ssROC is guaranteed to have lower variance than its supervised
counterpart, provided the PA score is not completely uninformative of the
gold-standard label. ssROC is also “safe” in the sense that it is equivalent
to supROC when the PA score is uninformative of the gold-standard label.

### Data and metrics for evaluation

We assessed the performance of ssROC using both simulated and real-world EHR
data from MGB.

#### Simulation study

We considered two settings that mimic scenarios with low and high PA accuracy.
In both settings, $Y$ is generated from a logistic regression model based on
$p$ covariates, ${\bf X}=(X_{1},\dots,X_{p})$. We generated ${\bf X}$ from a
multivariate normal distribution with each component of the mean vector
determined by a $\textrm{Bin}(3,0.3)$ random variable and a compound symmetric
covariance structure with a variance of 3 and covariance of 0.6. The two
settings are as follows:

  * •

Setting 1: PA with low accuracy. $Y$ is generated from the logistic regression
model $P(Y=1\mid{\bf
X})=\frac{1}{1+\exp\\{-(0.1X_{1}+0.1X_{2}+0.2X_{1}X_{2})\\}}$. The PA score is
generated using the linear predictor of the logistic regression model
excluding the interaction term so that
$S=\frac{1}{1+\exp\\{-(0.1X_{1}+0.1X_{2})\\}}$. In this setting, the AUC of
$S$ for $Y$ is 67.5.

  * •

Setting 2: PA with high accuracy. $Y$ is generated from the logistic
regression model $P(Y=1\mid{\bf
X})=\frac{1}{1+\exp\\{-(-4+X_{1}+X_{2}+0.5X_{3}+0.5X_{4})\\}}$. The PA score
is generated using the linear predictor of the logistic regression model so
that $S=\frac{1}{1+\exp\\{-(-4+X_{1}+X_{2}+0.5X_{3}+0.5X_{4})\\}}$. In this
setting, the AUC of $S$ for $Y$ is 95.3.

In both settings, we let the size of the unlabeled data be $N=10,000$ and
varied the labeled size as $n=100$ and $200$, resulting in a total of 4
simulation studies. We summarize results across 2,000 simulated datasets for
each study.

#### Real-world EHR data application

We validated our method using real-world EHR data from MGB, a Boston-based
healthcare system anchored by two tertiary care centers, Brigham and Women’s
Hospital and Massachusetts General Hospital. We evaluated the PAs for seven
phenotypes, including Breast Cancer (BrCa), Coronary Artery Disease (CAD),
Depression, Epilepsy, Hypertension (HTN), Type 1 Diabetes Mellitus (T1DM), and
Type 2 Diabetes Mellitus (T2DM) from the MGB Biobank. This biobank contains
linked EHR and genetic data from MGB between 1990 and 2015. We also evaluated
an algorithm for Pseudogout (pGout) from the Research Patient Data Registry of
MGB, which stores data on over 1 billion visits containing diagnoses,
medications, procedures, laboratory information, and clinical notes from 1991
to 2017.

The underlying full data for each phenotype consisted of patient records with
at least one phenotype-related PheCode in their record [47]. A subset of
patients was randomly sampled from the full data and sent for chart review.
The average size of the labeled and unlabeled data were 226 and 6,700,
respectively. For each phenotype, the PA was obtained by fitting PheNorm
without the random corruption denoising step. PheNorm is a recently proposed
weakly-supervised method based on normalizing silver-standard labels with
respect to patient healthcare utilization and is described in detail in
Supplementary Section 3 [48]. The silver-standard labels were taken as the
total numbers of (i) phenotype-related PheCodes and (ii) positive mentions of
the phenotype-related clinical concepts in patient records. Clinical concepts
were extracted from clinical notes with the Narrative Information Linear
Extraction and mapped to concept unique identifiers (CUIs) in the Unified
Medical Language System [49, 50]. The phenotypes represent varying levels of
PA accuracy, differences in labeled and unlabeled dataset sizes ($n$ and $N$,
respectively), and prevalence ($P$). A summary of the eight phenotypes is
presented in Table 1.

Phenotype | $n$ | $N$ | $P$ | PheCode | CUI |
---|---|---|---|---|---|---
Breast Cancer (BrCa) | 94 | 2002 | 0.78 | 174 | C0006142 |
Coronary Artery Disease (CAD) | 186 | 3793 | 0.37 | 411.4 | C0010054 |
Depression | 252 | 10,189 | 0.55 | 296.2 | C0011581 |
Epilepsy | 117 | 2225 | 0.48 | 345.1 | C0014544 |
Hypertension (HTN) | 390 | 19,853 | 0.79 | 401 | C0020538 |
pGout | 365 | 12,035 | 0.22 | 274 | C0263708 |
Type 1 Diabetes Mellitus (T1DM) | 128 | 2111 | 0.16 | 250.1 | C0011854 |
Type 2 Diabetes Mellitus (T2DM) | 280 | 3460 | 0.36 | 250.2 | C0011860 |

Table 1: Summary of the MGB Phenotypes. The labeled dataset size ($n)$, the
unlabeled dataset size ($N$), and the prevalence ($P$) of the eight phenotypes
as well as the main PheCode and concept unique identifier (CUI) used to train
PheNorm. The underlying full data for each phenotype included all participants
who passed the filter of $\geq 1$ PheCode for the phenotype of interest.

#### Benchmark method and reported metrics

We compared the PA evaluation results from ssROC and the benchmark, supROC,
using the real and simulated data. To improve the performance of the local
constant regression step, we transformed the PA scores by their empirical
cumulative distribution function [51]. We used the Gaussian kernel with the
bandwidth determined by the standard deviation of the standardized PA scores
divided by $n^{0.45}$ [38]. We obtained standard error estimates for the ROC
parameters using perturbation resampling with 500 replications and used
weights sampled from a scaled beta distribution given by
$4\cdot\text{Beta}(1/2,3/2)$ to improve finite-sample performance [52]. We
focused on logit-based confidence intervals (CIs) due to their improved
performance relative to standard Wald intervals [53]. To obtain the logit-
based intervals, we applied a logit transformation to the resampled estimates,
computed a Wald-based interval with the transformed estimates, and then
converted the upper and lower confidence bounds back to the original scale.

In simulations, the ground truth was the median of the ROC estimates using the
completely labeled data across the replicated datasets. We assessed percent
bias for both supROC and ssROC by computing the median of [(point estimate -
ground truth) / ground truth] $\times 100$ across the simulated datasets. We
also computed the empirical standard error (ESE) as the standard deviation of
the estimates across the replicated datasets and the asymptotic standard error
(ASE) as the median of the standard error estimates from the perturbation
resampling procedure across the replicated datasets. Using mean squared error
(MSE) as an aggregate measure of bias and variance, we evaluated the relative
efficiency (RE) as the ratio of the MSE of supROC to the MSE of ssROC. The
performance of our resampling procedure was assessed with the coverage
probability (CP) of the 95% CIs for both estimation procedures. In the real
data analysis, we present point estimates from both supROC and ssROC and the
RE defined as the ratio of the variance of supROC to ssROC. We evaluated the
performance of the PAs at an FPR of 10% and report the results for the AUC,
classification threshold (Threshold), TPR, PPV, and NPV for both the real and
simulated data.

## RESULTS

### Simulation

We present point estimates and percent biases for the ROC parameters at an FPR
of 10% for both simulation settings in Table 2. Both supROC and ssROC showed
minimal bias across all ROC parameters, labeled dataset sizes, and simulation
settings. The RE of ssROC compared to supROC is illustrated in Figure 4 with a
median RE of 1.2 for low PA accuracy (setting 1) and 1.5 for high PA accuracy
(setting 2). The higher gain in efficiency in setting 2 coincides with our
theoretical result and highlights the benefits of ssROC relative to supROC,
particularly when the PA is highly informative of the gold-standard label.

The most substantial improvements in the RE are observed for PPV. For example,
in setting 2 with $n=200$, the RE for PPV is 2.3. Practically, an RE exceeding
1 indicates that ssROC has a smaller variance than supROC (assuming the point
estimates from the two procedures are comparable). One may also view a
reduction in variance in terms of the sample size required for ssROC to
achieve the same variance as supROC. Again using PPV as an example, an RE of
2.3 suggests that ssROC with $n=87$ can achieve the same variance as supROC
with $n=200$. Supplementary Table 2 supports this observation, indicating that
the ASE of the PPV from ssROC with $n$ = 100 is lower than the ESE from supROC
with $n$ = 200 (2.0 vs. 2.2).

Supplementary Tables 1 and 2 present the ESE, ASE, and CP of the 95% CIs for
simulation settings 1 and 2, respectively. The proposed logit-based CI
consistently achieved coverage close to the nominal level for both methods
across all settings and parameters. The tables also illustrate the benefit of
the logit-based interval over the standard Wald interval. For example, in
setting 2 with $n$=100, the CPs of Wald CIs for AUC are 90.1% (supROC) and
89.5% (ssROC), while the CPs of the logit-based intervals are 96.9% and 93.5%.

|  | AUC | Threshold | TPR | PPV | NPV
---|---|---|---|---|---|---
$n$ | Method | Est. | %Bias | Est. | %Bias | Est. | %Bias | Est. | %Bias | Est. | %Bias
Setting 1
100 | supROC | 67.48 | -0.02 | 73.93 | -0.15 | 37.05 | 1.11 | 84.63 | 0.00 | 48.66 | 0.39
| ssROC | 67.55 | 0.14 | 73.84 | -0.14 | 37.18 | 1.13 | 84.59 | -0.03 | 48.84 | 0.62
200 | supROC | 67.45 | 0.04 | 73.87 | -0.01 | 36.90 | 0.70 | 84.33 | 0.03 | 48.70 | 0.24
| ssROC | 67.56 | 0.16 | 73.87 | -0.09 | 36.91 | 0.21 | 84.36 | 0.01 | 48.78 | 0.30
Setting 2
100 | supROC | 95.27 | 0.22 | 61.58 | 0.46 | 84.68 | 0.43 | 83.69 | -0.02 | 90.53 | 0.24
| ssROC | 95.10 | 0.02 | 61.55 | 0.36 | 84.36 | -0.51 | 83.79 | -0.14 | 90.37 | -0.21
200 | supROC | 95.30 | 0.12 | 61.48 | 0.25 | 85.08 | 0.47 | 83.85 | -0.04 | 90.77 | 0.25
| ssROC | 95.24 | 0.01 | 61.46 | 0.22 | 84.92 | 0.06 | 83.86 | -0.15 | 90.71 | 0.02

Table 2: Point estimates and percent bias for both simulation settings at a
FPR of 10%. For both settings, the size of the unlabeled was $N=10,000$.
Figure 4: Relative efficiency for both simulation settings at a FPR of 10%.
Relative efficiency is defined as the mean squared error of supROC compared to
the mean squared error of ssROC. For both settings, the size of the unlabeled
was $N=10,000$.

### Analysis of eight PAs from MGB

Figure 5 displays the point estimates for the eight phenotypes from MGB at a
FPR of 10%. The top row contains the PAs with lower accuracy. Generally, the
point estimates from ssROC are similar to those from supROC. However, there
are some differences in the classification threshold estimates for Depression,
HTN, and CAD, which leads to some discrepancies in the corresponding point
estimates for the TPR, PPV, and NPV. As supROC is only evaluated at the unique
PA scores in the labeled dataset, the threshold estimate can be unstable at
some FPRs. In contrast, ssROC evaluates the ROC parameters across a broader
range of PA scores in the larger unlabeled dataset and results in a more
stable estimation. Despite these differences, it is important to note that the
95% CIs for all estimators consistently overlapped, as shown in Supplementary
Figure 1.

Figure 6 shows the median RE of supROC to ssROC across the eight phenotypes at
a FPR of 10%. The median RE gain across phenotypes and metrics was 1.7,
implying that ssROC can achieve the same precision as supROC, but with roughly
60% of the amount of data on average. Consistent with the simulation findings,
the most significant RE gain was observed in PPV followed by the
classification threshold, with median RE both near 2.1. The lowest gain was in
AUC, with a median RE of 1.2. These results suggest that ssROC can achieve
equivalent precision to supROC with respect to PPV with approximately 52% of
the labeled data and maintain equivalent precision for AUC with 15% less
labeled data. Additionally, Depression and HTN exhibited very high RE for the
classification threshold (Supplementary Figure 2). This behavior further
emphasizes the stability of ssROC relative to supROC in real data. For
instance, Supplementary Figure 3 illustrates that the empirical distribution
of the resampled estimates from supROC is skewed for Depression and multimodal
for HTN due to the aforementioned difficulties in estimating the
classification threshold with the small amount of labeled data.

Consistent with our simulation and theoretical results, we also observe that
RE is linked to PA accuracy. For example, a phenotype with high PA accuracy,
such as BrCa (AUC = 94.4), has a median RE of 2.2 while pGout, an episodic
disease without a specific billing code, has the lowest PA accuracy (AUC =
61.4) and a median RE near 1. These findings underscore the advantages of our
proposed ssROC method compared to its supervised counterpart: (i) ssROC
increases precision when the PA score is predictive and (ii) performs on par
with the supROC when the PA score is uninformative.

Figure 5: Point estimates for the 8 phenotypes from MGB at a FPR of 10%. Point
estimates of the ROC parameters for both ssROC and supROC. Figure 6: Relative
efficiency for the 8 phenotypes from MGB at a FPR of 10%. The median relative
efficiency (RE) of supROC versus ssROC across the 8 phenotypes. RE is defined
as the ratio of the variance of the estimates from supROC and the variance of
the estimates from ssROC.

## DISCUSSION

Although high-throughput phenotyping is the backbone of EHR-based research,
there is a paucity of methods for reliably evaluating the predictive
performance of a PA with limited labeled data. The proposed ssROC method fills
this gap. ssROC is a simple two-step estimation procedure that leverages large
volumes of unlabeled data by imputing missing gold-standard labels with a
nonparametric recalibration of a PA score. Unlike existing procedures for PA
evaluation in the informatics literature, ssROC eliminates the requirement
that the PA be properly calibrated to yield unbiased estimation of the ROC
parameters and may be utilized for ML-based PAs [25, 37]. While we focus
specifically on weakly-supervised PAs in our theoretical analysis and real
data example, given their increasing popularity and ability to automate PA
estimation, ssROC can also be used to evaluate rule-based or other ML-based
PAs. Moreover, by harnessing unlabeled data, ssROC yields substantially less
variable estimates than supROC in both simulated and real data. Practically,
this translates into a significant reduction in the amount of chart review
required to obtain a precise understanding of PA performance.

While our work is a first step toward streamlining PA evaluation, there are
several avenues for future research. First, ssROC assumes that the labeled
examples have been randomly sampled from the underlying full data. More
effective sampling strategies, such as targeted or stratified sampling, have
the potential to further improve upon the efficiency of ssROC, but require
future research [54]. Secondly, the nonparametric recalibration step demands a
sufficient amount of labeled data for the local constant regression to be well
estimated. Our future work will develop a flexible parametric recalibration
procedure that accommodates smaller labeled data sizes. Thirdly, ssROC can
also be extended for model comparisons and evaluation of fairness metrics,
which are urgently needed given the increasing recognition of algorithmic
unfairness in informatics applications. The recalibration step would need to
be augmented in both settings to utilize additional information in multiple PA
scores or protected attributes, respectively. This augmentation could lead to
a potentially more efficient procedure as ssROC currently only uses
information from one PA score for imputation. Lastly, our results demonstrate
the ability of ssROC to provide accurate phenotyping for eight phenotypes with
variable prevalence, labeled and unlabeled dataset sizes, and PA performance
within one health system. Further work is needed to understand the performance
of our method across a diverse range of phenotypes and to extend our approach
to accommodate federated analyses across multiple healthcare systems.

## CONCLUSION

In this paper, we introduced a semi-supervised approach, called ssROC, that
leverages a large volume of unlabeled data together with a small subset of
gold-standard labeled to precisely estimate the ROC parameters of PAs. PA
development involves two key steps: (i) algorithm estimation and (ii)
algorithm evaluation. While a considerable amount of effort has been placed on
algorithm estimation, ssROC fills the current gap in robust methodology for
predictive performance evaluation. Additionally, ssROC is simple to implement
and is available in open-source R software to encourage use in practice. When
used in conjunction with weakly-supervised PAs, ssROC demonstrates the
potential to facilitate the reliable and streamlined phenotyping that is
necessary for a wide variety of translational EHR applications.

## FUNDING

The project was supported by the Natural Sciences and Engineering Research
Council of Canada grant (RGPIN-2021-03734), the University of Toronto
Connaught New Researcher Award, and the University of Toronto Seed Funding for
Methodologists Grant (to JeG).

## AUTHOR CONTRIBUTIONS

JeG conceived and designed the study. JeG, JiG, CB, and CH conducted
simulation analyses. CB and CH conducted real data analyses. JeG, JiG, CB, and
KZ analyzed and interpreted the results. JeG and JiG drafted and revised the
manuscript. All authors reviewed and approved the final manuscript.

## CONFLICT OF INTEREST STATEMENT

The authors have no conflicts of interest to declare.

## DATA AND CODE AVAILABILITY

Our proposed method is implemented as an R software package, ssROC, which is
available at https://github.com/jlgrons/ssROC.

## References

  * [1] J Michael McGinnis et al.  “Clinical data as the basic staple of health learning: Creating and protecting a public good: Workshop summary”  National Academies Press, 2011
  * [2] Kenneth S Boockvar et al.  “Electronic health records and adverse drug events after patient transfer”  In _Quality and Safety in Health Care_ 19.5 BMJ Publishing Group Ltd, 2010, pp. e16–e16
  * [3] Fina Kurreeman et al.  “Genetic basis of autoantibody positive and negative rheumatoid arthritis risk in a multi-ethnic cohort derived from electronic health records”  In _The American Journal of Human Genetics_ 88.1 Elsevier, 2011, pp. 57–69
  * [4] Katherine P Liao et al.  “Associations of autoantibodies, autoimmune risk alleles, and clinical diagnoses from the electronic medical records in rheumatoid arthritis cases and non–rheumatoid arthritis controls”  In _Arthritis & Rheumatism_ 65.3 Wiley Online Library, 2013, pp. 571–581
  * [5] Chia-Yen Chen et al.  “Genetic validation of bipolar disorder identified by automated phenotyping using electronic health records”  In _Translational psychiatry_ 8.1 Nature Publishing Group, 2018, pp. 1–8
  * [6] Ruowang Li et al.  “Electronic health records and polygenic risk scores for predicting disease risk”  In _Nature Reviews Genetics_ Nature Publishing Group, 2020, pp. 1–10
  * [7] Gabriel A Brat et al.  “International electronic health record-derived COVID-19 clinical course profiles: the 4CE consortium”  In _NPJ digital medicine_ 3.1 Nature Publishing Group, 2020, pp. 1–9
  * [8] Lisa Bastarache  “Using phecodes for research with the electronic health record: from PheWAS to PheRS”  In _Annual review of biomedical data science_ 4 Annual Reviews, 2021, pp. 1–19
  * [9] Daniel Prieto-Alhambra et al.  “Unraveling COVID-19: a large-scale characterization of 4.5 million COVID-19 cases using CHARYBDIS”  In _Research square_ American Journal Experts, 2021
  * [10] Katharine E Henry et al.  “Factors driving provider adoption of the TREWS machine learning-based early warning system and its effects on sepsis treatment timing”  In _Nature medicine_ 28.7 Nature Publishing Group US New York, 2022, pp. 1447–1454
  * [11] Chaitanya Shivade et al.  “A review of approaches to identifying patient phenotype cohorts using electronic health records”  In _Journal of the American Medical Informatics Association_ 21.2 BMJ Publishing Group, 2014, pp. 221–230
  * [12] Juan M Banda et al.  “Advances in electronic phenotyping: from rule-based definitions to machine learning models”  In _Annual review of biomedical data science_ 1 NIH Public Access, 2018, pp. 53
  * [13] Hadeel Alzoubi et al.  “A review of automatic phenotyping approaches using electronic health records”  In _Electronics_ 8.11 MDPI, 2019, pp. 1235
  * [14] Siyue Yang et al.  “Machine learning approaches for electronic health records phenotyping: a methodical review”  In _Journal of the American Medical Informatics Association_ 30.2 Oxford University Press, 2023, pp. 367–381
  * [15] Yichi Zhang et al.  “High-throughput phenotyping with electronic medical record data using a common semi-supervised approach (PheCAP)”  In _Nature Protocols_ 14.12 Nature Publishing Group, 2019, pp. 3426–3444
  * [16] Shawn Murphy et al.  “Instrumenting the health care enterprise for discovery research in the genomic era”  In _Genome research_ 19.9 Cold Spring Harbor Lab, 2009, pp. 1675–1681
  * [17] Victor Castro et al.  “Identification of subjects with polycystic ovary syndrome using electronic health records”  In _Reproductive Biology and Endocrinology_ 13.1 Springer, 2015, pp. 116
  * [18] Pedro L Teixeira et al.  “Evaluating electronic health record data sources and algorithmic approaches to identify hypertensive individuals”  In _Journal of the American Medical Informatics Association_ 24.1, 2017, pp. 162–171
  * [19] Alon Geva et al.  “A computable phenotype improves cohort ascertainment in a pediatric pulmonary hypertension registry”  In _The Journal of pediatrics_ 188 Elsevier, 2017, pp. 224–231
  * [20] Christopher Meaney et al.  “Using Biomedical Text as Data and Representation Learning for Identifying Patients with an Osteoarthritis Phenotype in the Electronic Medical Record”  In _International Journal of Population Data Science_ 3.4, 2018
  * [21] Sebastian Gehrmann et al.  “Comparing deep learning and concept extraction based methods for patient phenotyping from clinical narratives”  In _PLOS ONE_ 13, 2018, pp. e0192360
  * [22] Katherine P Liao et al.  “High-throughput multimodal automated phenotyping (MAP) with application to PheWAS”  In _Journal of the American Medical Informatics Association_ 26.11 Oxford University Press, 2019, pp. 1255–1262
  * [23] Vijay S. Nori et al.  “Deep neural network models for identifying incident dementia using claims and EHR datasets”  In _PLOS ONE_ 15, 2020, pp. e0236400
  * [24] Yizhao Ni et al.  “Automated detection of substance use information from electronic health records for a pediatric population”  In _Journal of the American Medical Informatics Association_ 28, 2021, pp. 2116–2127
  * [25] Joel N Swerdel, George Hripcsak and Patrick B Ryan  “PheValuator: development and evaluation of a phenotype algorithm evaluator”  In _Journal of biomedical informatics_ 97 Elsevier, 2019, pp. 103258
  * [26] Christian Chartier, Lisa Gfrerer and William G. Austen  “ChartSweep: A HIPAA-compliant Tool to Automate Chart Review for Plastic Surgery Research”  In _Plastic and Reconstructive Surgery - Global Open_ 9, 2021, pp. e3633
  * [27] Sheng Yu et al.  “Toward high-throughput phenotyping: unbiased automated feature extraction and selection from knowledge sources”  In _Journal of the American Medical Informatics Association_ 22.5 The Oxford University Press, 2015, pp. 993–1000
  * [28] Sheng Yu et al.  “Surrogate-assisted feature extraction for high-throughput phenotyping”  In _Journal of the American Medical Informatics Association_ 24.e1 Oxford University Press, 2017, pp. e143–e149
  * [29] Isabelle-Emmanuella Nogues et al.  “Weakly Semi-supervised phenotyping using Electronic Health records”  In _Journal of Biomedical Informatics_ 134 Elsevier, 2022, pp. 104175
  * [30] Adam Wright, Elizabeth S Chen and Francine L Maloney  “An automated technique for identifying associations between medications, laboratory results and problems”  In _Journal of biomedical informatics_ 43.6 Elsevier, 2010, pp. 891–901
  * [31] Adam Wright et al.  “A method and knowledge base for automated inference of patient problems from structured data in an electronic medical record”  In _Journal of the American Medical Informatics Association_ 18.6 BMJ Group BMA House, Tavistock Square, London, WC1H 9JR, 2011, pp. 859–867
  * [32] Vibhu Agarwal et al.  “Learning statistical models of phenotypes using noisy labeled training data”  In _Journal of the American Medical Informatics Association_ 23.6 Oxford University Press, 2016, pp. 1166–1173
  * [33] Juan M Banda et al.  “Electronic phenotyping with APHRODITE and the Observational Health Sciences and Informatics (OHDSI) data network”  In _AMIA Summits on Translational Science Proceedings_ 2017 American Medical Informatics Association, 2017, pp. 48
  * [34] Jing Huang et al.  “PIE: A prior knowledge guided integrated likelihood estimation method for bias reduction in association studies using electronic health records data”  In _Journal of the American Medical Informatics Association_ 25.3 Oxford University Press, 2018, pp. 345–352
  * [35] Jiayi Tong et al.  “An augmented estimation procedure for EHR-based association studies accounting for differential misclassification”  In _Journal of the American Medical Informatics Association_ 27.2 Oxford University Press, 2020, pp. 244–253
  * [36] Ziyan Yin et al.  “A cost-effective chart review sampling design to account for phenotyping error in electronic health records (EHR) data”  In _Journal of the American Medical Informatics Association_ 29.1 Oxford University Press, 2022, pp. 52–61
  * [37] Joel N Swerdel et al.  “PheValuator 2.0: Methodological improvements for the PheValuator approach to semi-automated phenotype algorithm evaluation”  In _Journal of Biomedical Informatics_ 135 Elsevier, 2022, pp. 104177
  * [38] Jessica L Gronsbell and Tianxi Cai  “Semi-supervised approaches to efficient evaluation of model prediction performance”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 80.3 Wiley Online Library, 2018, pp. 579–594
  * [39] Jessica Gronsbell et al.  “Efficient estimation and evaluation of prediction rules in semi-supervised settings under stratified sampling”  In _arXiv preprint arXiv:2010.09443_ , 2020
  * [40] Ben Van Calster et al.  “Calibration: the Achilles heel of predictive analytics”  In _BMC medicine_ 17.1 Springer, 2019, pp. 1–7
  * [41] Yingxiang Huang et al.  “A tutorial on calibration measurements and calibration models for clinical prediction models”  In _Journal of the American Medical Informatics Association_ 27.4 Oxford University Press, 2020, pp. 621–633
  * [42] Margaret Sullivan Pepe  “The statistical evaluation of medical tests for classification and prediction”  Oxford University Press, 2003
  * [43] Tianxi Cai and Chaya S Moskowitz  “Semi-parametric estimation of the binormal ROC curve for a continuous diagnostic test”  In _Biostatistics_ 5.4 Oxford University Press, 2004, pp. 573–586
  * [44] Zhezhen Jin, Zhiliang Ying and L.J. Wei  “A Simple Resampling Method by Perturbing the Minimand”  In _Biometrika_ 88.2, 2001, pp. 381–390
  * [45] Patrice Bertail, Stéphan Clémençcon and Nicolas Vayatis  “On bootstrapping the ROC curve”  In _Advances in Neural Information Processing Systems_ 21, 2008
  * [46] Whitney Newey, Fushing Hsieh and James Robins  “Undersmoothing and Bias Corrected Functional Estimation”, 1998  URL: https://EconPapers.repec.org/RePEc:mit:worpap:98-17
  * [47] Joshua C Denny et al.  “Systematic comparison of phenome-wide association study of electronic medical record data and genome-wide association study data”  In _Nature biotechnology_ 31.12 Nature Publishing Group, 2013, pp. 1102–1111
  * [48] Sheng Yu et al.  “Enabling phenotypic big data with PheNorm”  In _Journal of the American Medical Informatics Association_ 25.1 Oxford University Press, 2018, pp. 54–60
  * [49] Donald AB Lindberg, Betsy L Humphreys and Alexa T McCray  “The unified medical language system”  In _Yearbook of medical informatics_ 2.01 Georg Thieme Verlag KG, 1993, pp. 41–51
  * [50] Sheng Yu, Tianrun Cai and Tianxi Cai  “NILE: fast natural language processing for electronic health records”  In _arXiv preprint arXiv:1311.6063_ , 2013
  * [51] Matthew P Wand, James Stephen Marron and David Ruppert  “Transformations in density estimation”  In _Journal of the American Statistical Association_ 86.414 Taylor & Francis, 1991, pp. 343–353
  * [52] Jennifer A Sinnott and Tianxi Cai  “Inference for survival prediction under the regularized Cox model”  In _Biostatistics_ 17.4 Oxford University Press, 2016, pp. 692–707
  * [53] Alan Agresti  “Categorical data analysis”  John Wiley & Sons, 2012
  * [54] W Katherine Tan and Patrick J Heagerty  “Surrogate-guided sampling designs for classification of rare outcomes from electronic medical records data”  In _Biostatistics_ 23.2 Oxford University Press, 2022, pp. 345–361


# Scaling Graph-based Deep Learning models to larger networks

Miquel Ferriol-Galmés1, José Suárez-Varela1, Krzysztof Rusek1,2, Pere Barlet-
Ros1, Albert Cabellos-Aparicio1 1Barcelona Neural Networking center,
Universitat Politècnica de Catalunya, Spain 2 AGH University of Science and
Technology, Krakow, Poland

###### Abstract

Graph Neural Networks (GNN) have shown a strong potential to be integrated
into commercial products for network control and management. Early works using
GNN have demonstrated an unprecedented capability to learn from different
network characteristics that are fundamentally represented as graphs, such as
the topology, the routing configuration, or the traffic that flows along a
series of nodes in the network. In contrast to previous solutions based on
Machine Learning (ML), GNN enables to produce accurate predictions even in
other networks unseen during the training phase. Nowadays, GNN is a hot topic
in the Machine Learning field and, as such, we are witnessing great efforts to
leverage its potential in many different fields (e.g., chemistry, physics,
social networks). In this context, the Graph Neural Networking challenge 2021
brings a practical limitation of existing GNN-based solutions for networking:
the lack of generalization to larger networks. This paper approaches the
scalability problem by presenting a GNN-based solution that can effectively
scale to larger networks including higher link capacities and aggregated
traffic on links.

###### Index Terms:

Network Modeling, Graph Neural Networks.

##  I Introduction

Graph Neural Networks (GNN) have produced groundbreaking applications in many
fields where data is fundamentally structured as graphs (e.g., chemistry,
physics, biology, recommender systems). In the field of computer networks,
this new type of neural networks is being rapidly adopted for a wide variety
of use cases [1], particularly for those involving complex inter-dependencies
between different network elements (e.g., performance modeling, routing
optimization, resource allocation in wireless networks). In the context of
network modeling, unlike previous solutions based on Machine Learning (ML),
GNN enables to produce accurate predictions even in networks unseen during the
training phase. Nowadays, GNN is a hot topic in the ML field and, as such, we
are witnessing great efforts to leverage its potential in many different
fields where data is fundamentally represented as graphs (e.g., chemistry,
physics, social networks). The Graph Neural Networking challenge [2] is an
annual competition that brings fundamental challenges on the application of
GNN to networking.

The 2021 edition of the Graph Neural Networking challenge [2] brings a
fundamental limitation of existing GNNs: their lack of generalization
capability to larger graphs. In order to achieve production-ready GNN-based
solutions, we need models that can be trained in network testbeds of limited
size (e.g., at the vendor’s networking lab), and then be directly ready to
operate with guarantees in real customer networks, which are often much larger
in number of nodes. In this challenge, participants are asked to design GNN-
based models that can be trained on small network scenarios (up to 50 nodes),
and after that scale successfully to larger networks not seen before, up to
300 nodes.

This paper presents a possible solution for the Graph Neural Networking
challenge 2021, which overcomes the scalability limitations of current GNN
models applied to networking. The results show that the proposed model is able
to accurately generalize from small-scale networks seen during training (up to
50 nodes) to samples of networks considerably larger (up to 300 nodes) unseen
in advance by the model.

##  II Problem statement

The goal of this challenge is to create a scalable Network Digital Twin (i.e.,
a network model) based on neural networks, which can accurately estimate QoS
performance metrics given a network state snapshot. More in detail, solutions
must predict the resulting source-destination mean per-packet delay given:
$(i)$ a network topology, $(ii)$ a source-destination traffic matrix, and
$(iii)$ a routing configuration (see Figure 1).

Figure 1: Scheme of the neural network-based solution requested in the Graph
Neural Networking challenge [2]

##  III Background: GNNs

Graphs are used to represent relational information. Particularly, a graph
$G\in\\{V,E\\}$ comprises a set of objects $V$ (i.e., vertices) and some
relations between them $E$ (i.e., edges).

GNN [3] is a family of NNs especially designed to work with graph-structured
data. These models dynamically build their internal NN architecture based on
the input graph. For this, they use a modular NN structure that represents
explicitly the elements and connections of the graph. As a result, they
support graphs of variable size and structure, and their graph processing
mechanism is equivariant to node and edge permutation, which eventually endows
them with strong generalization capabilities over graphs – also known as
strong relational inductive bias [4].

Despite GNN covers a broad family of neural networks with different
architectural variants (e.g., [3, 5, 6]), most of them share the basic
principle of an iterative message-passing phase, where the different elements
of the graph exchange information according to their connections, and a final
readout phase uses the information encoded in graph elements to produce the
final output(s). We refer the reader to [3, 7, 4] for a more comprehensive
background on GNNs.

In the context of computer networks, standard GNN architectures, as the one
described above, are not applicable, as networks introduce graphs with
heterogeneous elements and complex circular dependencies between them. This
requires devising a more complex GNN architecture that can adapt to the
intricacies of computer networks.

##  IV GNN-Based Solution

This section describes a novel GNN-based solution tailored to accurately model
the behavior of real network infrastructures. The model implements a novel
three-stage message passing algorithm that explicitly defines some key
elements for network modeling (e.g. forwarding devices, queues, paths), and
offers support for a wide variety of features introduced in modern networking
trends (e.g., complex QoS-aware queuing policies, overlay routing).

Figure 2: Back-box scheme of the proposed GNN-based model.

Figure 2 shows a black-box representation of the proposed GNN-based network
model. The input of this model is a network state sample, defined by: a
network topology, a set of traffic models (flow-level), a routing scheme
(flow-level), and a queuing configuration (interface-level). As output, this
model produces estimates of relevant performance metrics at a flow-level
granularity (e.g., delay, jitter, losses). Note that, beyond the solution
requested in the Graph Neural Networking challenge (Fig. 1) the proposed model
also supports different queueing policies (e.g., Weighted Fair Queuing,
Deficit Round Robin), while in the challenge all forwarding devices implement
a FIFO policy.

###  IV-A Model Description

The proposed GNN-Based solution pursues two main objectives: $(1)$ Finding a
good representation for the network components supported by the model (e.g.,
traffic intensities, routing), and $(2)$ Exploit scale-independent features of
networks, in order to achieve good generalization to larger networks than
those seen during training.

1) Representing network components and their relationships:

First, let us define a network as a set of links
$L=\\{l_{i}:i\in(1,...,n_{l})\\}$, a set of queues on
$Q=\\{q_{i}:i\in(1,...,n_{q})\\}$, and a set of source-destination flows
$F=\\{f_{i}:i\in(1,...,n_{f})\\}$. According to the routing configuration,
flows follow a source-destination path. Hence, we define flows as sequences
with tuples of the queues and links they traverse
$f_{i}=\\{(q_{F_{q}(f_{i},0)},l_{F_{l}(f_{i},0)}),...,(q_{F_{q}(f_{i},|f_{i}|)},$
$l_{F_{l}(f_{i},|f_{i}|)})\\}$, where $F_{q}(f_{i},j)$ and $F_{l}(f_{i},j)$
respectively return the index of the $j$-th queue or link along the path of
flow $f_{i}$. Let us also define $Q_{f}(q_{i})$ as a function that returns all
the flows passing through queue $q_{i}$, and $L_{q}(l_{i})$ as a function that
returns the queues injecting traffic into link $l_{i}$ – i.e., the queues at
the output port to which the link is connected.

Figure 3: Schematic representation of the network model implemented.

Following the previous notation, the GNN-based model considers an input graph
with three main components: $(i)$ the physical links $L$ that shape the
network topology, $(ii)$ the queues $Q$ at each output port of network
devices, and $(iii)$ the active flows $F$ in the network, which follow some
specific src-dst paths (i.e., sequences of queues and links), and whose
traffic is generated from a given traffic model. Figure 3 shows a schematic
representation of the network model internally considered, which is derived
from the several mechanisms that affect performance in real networks. From
this model, we can extract three basic principles:

  1. (i)

The state of flows (e.g., throughput, losses) is affected by the state of the
queues and links they traverse (e.g., queue/link utilization).

  2. (ii)

The state of queues (e.g., occupation) depends on the state of the flows
passing through them (e.g., traffic model).

  3. (iii)

The state of links (e.g., utilization) depends on the states of the queues at
the output port of the link.

Formally, these principles can be formulated as follows:

$\displaystyle h_{f_{k}}$
$\displaystyle=g_{f}(h_{q_{k(0)}},h_{l_{k(0)}},...,h_{q_{k(|f_{k}|)}},h_{l_{k(|f_{k}|)}})$
(1) $\displaystyle h_{q_{i}}$
$\displaystyle=g_{q}(h_{p_{1}},...,h_{p_{m}}),\quad q_{i}\in p_{k},k=1,...,j$
(2) $\displaystyle h_{l_{j}}$
$\displaystyle=g_{l}(h_{q_{1}},...,h_{q_{m}}),\quad q_{m}\in L_{q}(l_{j})$ (3)

Where $g_{f}$, $g_{q}$ and $g_{l}$ are some unknown functions, and $h_{f}$,
$h_{q}$ and $h_{l}$ are latent variables that encode information about the
state of flows $F$, queues $Q$, and links $L$ respectively. Note that these
principles define a circular dependency between the three network components
($F$, $Q$, and $L$) that must be solved to find latent representations
satisfying the equations above.

Based on the previous network modeling principles, we define the architecture
of the model (see Algorithm 1). Our GNN-based model implements a custom three-
stage message-passing algorithm that combines the states of flows, queues and
links according to Equations (1)-(3), thus aiming to resolve the circular
dependencies defined in such functions. First, the hidden states $h_{l}$,
$h_{q}$, and $h_{f}$ – represented as n-element vectors – are initialized with
some features (lines 3-5), denoted respectively by $x_{l_{i}}$, $x_{q_{j}}$
and $x_{f_{k}}$. In our case, we set the initial features of links ($x_{l}$)
as: $(i)$ the link capacity ($C_{i}$), and $(ii)$ the scheduling policy at the
output port of the link (FIFO, SP, WFQ, or DRR [8]), using one-hot encoding.
For the initial features of queues ($x_{q}$) we include: $(i)$ the buffer
size, $(ii)$ the priority level (one-hot encoding), and $(iii)$ the weight
(only for WFQ and DRR). Lastly, the initial flow features ($x_{f}$) are a
descriptor of the traffic model used in the flow ($T_{i}$). Once the states
are initialized, the message-passing phase is iteratively executed $T$ times
(loop from line 6), where $T$ is a configurable parameter. Each message-
passing iteration is in turn divided in three stages, that respectively
represent the message passing and update of the hidden states of flows $h_{f}$
(lines 7-11), queues $h_{q}$ (lines 12-15), and links $h_{l}$ (lines 16-19).

Algorithm 1 Internal architecture of the proposed model

1:$F$, $Q$, $L$, $x_{f}$, $x_{q}$, $x_{l}$

2:$h^{T}_{q}$, $h^{T}_{l}$, $h^{T}_{f}$, $\hat{y_{f}}$, $\hat{y_{q}}$,
$\hat{y_{l}}$

3:for each $l\in L$ do $h^{0}_{l}\leftarrow[x_{l},0...0]$

4:for each $q\in Q$ do $h^{0}_{q}\leftarrow[x_{q},0...0]$

5:for each $f\in F$ do $h^{0}_{f}\leftarrow[x_{f},0...0]$

6:for t = 0 to T-1 do $\triangleright$ Message Passing Phase

7: for each $f\in F$ do $\triangleright$ Message Passing on Flows

8: for each $(q,l)\in f$ do

9: $h^{t}_{f}\leftarrow FRNN(h^{t}_{f},[h^{t}_{q},h^{t}_{l}])$
$\triangleright$ Flow: Aggr. and Update

10: $\widetilde{m}^{t+1}_{f,q}\leftarrow h^{t}_{f}$ $\triangleright$ Flow:
Message Generation

11: $h^{t+1}_{f}\leftarrow h^{t}_{f}$

12: for each $q\in Q$ do $\triangleright$ Message Passing on Queues

13: $M_{q}^{t+1}\leftarrow\sum_{f\in Q_{f}(q)}\widetilde{m}^{t+1}_{f,q}$
$\triangleright$ Queue: Aggregation

14: $h^{t+1}_{q}\leftarrow U_{q}(h^{t}_{q},M_{q}^{t+1})$ $\triangleright$
Queue: Update

15: $\widetilde{m}^{t+1}_{q}\leftarrow h^{t+1}_{q}$ $\triangleright$ Queue:
Message Generation

16: for each $l\in L$ do $\triangleright$ Message Passing on Links

17: for each $q\in L_{q}(l)$ do

18: ${h}^{t}_{l}\leftarrow LRNN(h^{t}_{l},\widetilde{m}^{t+1}_{q})$
$\triangleright$ Link: Aggr. and Update

19: $h^{t+1}_{l}\leftarrow{h}^{t}_{l}$

20:$\hat{y_{f}}\leftarrow R_{f}(h^{T}_{f})$ $\triangleright$ Readout phase

21:$\hat{y_{q}}\leftarrow R_{q}(h^{T}_{q})$

Finally, functions $R_{f}$ (line 20) and $R_{q}$ (line 21) represent
independent readout functions that can be respectively applied to the hidden
states of flows $h_{f}$ or queues $h_{q}$.

2) Scaling to larger networks: scale-independent features

GNNs have shown an unprecedented capability to generalize over graph-
structured data [4, 9]. In the context of generalizing to larger graphs, it is
well known that these models keep good generalization capabilities as long as
the spectral properties of graphs are similar to those seen during training
[10]. In the case of the proposed model, its message-passing algorithm can
analogously generalize to graphs with similar structures to those seen during
the training phase – e.g., graphs with a similar number of queues at output
ports, or a similar number of flows aggregated in queues. In this vein,
generating a representative dataset in small networks, covering a wide range
of graph structures, does not imply any practical limitation to then achieve
good generalization properties to larger networks. It can be done by simply
adding a broad combination of realistic network samples (e.g., with a wide
variety of traffic models, routing schemes, queuing policies) – as in the
process described later in section IV-B.

However, from a practical standpoint, scaling to larger networks often entails
a broader definition beyond the topology size and structure. In particular,
there are two main properties we can observe as networks become larger: $(i)$
higher link capacities (as there is more aggregated traffic in the core links
of the network), and $(ii)$ longer paths (as the network diameter becomes
larger). This requires devising mechanisms to effectively scale on these two
features.

Scaling to larger link capacities: If we observe the internal architecture of
the model (Algorithm 1), we can find that the link capacity $C$ is only
represented as an initial feature of links’ hidden states $x_{l_{i}}$. The
fact that $C$ is encoded as a numerical feature in the model introduces
inherent limitations to scale to larger capacity values. Indeed, scaling to
out-of-distribution numerical values is widely recognized as a generalized
limiting factor among all neural networks [11, 12]. Thus, our approach is to
exploit particularities from the network domain to find scale-independent
representations that can define link capacities and how they relate to other
link-level features that impact on performance (e.g., the aggregated traffic
in the link), as the final goal of this work is to accurately estimate
performance metrics (e.g., delay, jitter, losses). Inspired by traditional
queuing theory (QT) methods, we aim to encode the relative ratio between the
arrival rates on links (based on the traffic aggregated in the link), and the
service times (based on the link capacity), thus enabling the possibility to
infer the output performance metrics of our model from scale-independent
values. As a result, we define link capacities ($Cap_{link}$) as the product
of a virtual reference link capacity ($C_{ref}$) and a scale factor ($S_{f}$):

$Cap_{link}=C_{ref}*S_{f}$ (4)

This representation enables to define arbitrary combinations of scale factors
and reference link capacities to define the actual capacity of links in
networks. Hence, we introduce the capacity feature ($C_{i}$) as a 2-element
vector defined as $C_{i}$$=$$[C_{ref},S_{f}]$, which is included in the
initial feature vector of links ($x_{l}$). Note that this feature will
eventually be encoded in the hidden states of links ($h_{l}$). In the internal
architecture of the GNN-based model (Algorithm 1), this factor will mainly
affect the update functions of flows and links (lines 9 and 18), as they are
the only ones that process directly the hidden states of links ($h_{l}$). As a
result, the RNNs approximating these update functions can potentially learn to
make accurate estimates on any combination of $C_{ref}$ and $S_{f}$ as long as
these two features are within the range of values observed independently for
each of them during the training phase (i.e.,
$S_{f}\in[s_{f_{min}},s_{f_{max}}]$ and
$C_{ref}\in[C_{{ref}_{min}},C_{{ref}_{max}}]$). Thus, we exploit this property
to devise a custom data augmentation method, where we take samples from small
networks with limited link capacities and generate different combinations of
$C_{ref}$ and $S_{factor}$ that enable us to scale accurately to considerably
larger capacities. Note that in this process, the numerical values seen by the
model ($C_{ref}$ and $S_{factor}$) are kept in the same ranges both in the
training on small networks and the posterior inference on larger networks,
thus overcoming the practical limitation of out-of-distribution predictions
[11, 12]. More details about the proposed data augmentation process are given
in Sec. IV-C.

The previous mechanism enables to keep scale-independent features along the
message-passing phase of our model (loop lines 6-19 in Algorithm 1), while it
is still needed to extend the scale independence to the output layer of the
model. Particularly, in this paper, we use the model to predict the flows’
delays. Note that the distribution of these parameters can also vary for flows
traversing links with higher capacities, thus leading again to out-of-
distribution values. Based on the fundamentals of QT, we overcome this
potential limitation by inferring delays indirectly from the occupation of
queues in the network $O_{q_{i}}$$\in$$[0,1]$, using the
$\hat{y_{q}}$$=$$R_{q}(h_{q})$ function of Algorithm 1. Then, we infer the
flow delay as a linear combination of the waiting times in queues (inferred
from $O_{q_{i}}$) and the transmission times of the links the flow traverses.
Note that a potential advantage with respect to traditional QT models is that
the queue occupation estimates produced can be more accurate, especially for
complex traffic models resembling real-world traffic.

Scaling to longer paths: In the internal architecture of the proposed GNN, the
path length only affects to the RNN function of line 9 (Algorithm 1), which
collects the state of queues ($h_{q}$) and links ($h_{l}$) to update flows’
states ($h_{f}$). The main limitation here is that this RNN can typically see
during training shorter link-queue sequences than those it can find then in
larger networks, which can potentially have longer paths. As a result, we
define $L_{max}$ as a configurable parameter of our model that defines the
maximum sequence length supported by this RNN. Then, we split flows exceeding
$L_{max}$ into different queue-link sequences that are independently digested
by the RNN. To keep the state along the whole flow, in case it is divided into
more than one sequence, we initialize the initial state of the RNN with the
output resulting from the previous sequence.

###  IV-B Simulation Setup

To train, validate and test the model we use as ground truth a packet-level
network simulator (OMNeT++ v5.5.1 [13]). Each sample is labeled with network
performance metrics obtained by the simulator: per-source-destination
performance measurements (mean per-packet delay, jitter and loss), and port
statistics (e.g., queue utilization, size). To generate these datasets, for
each sample we randomly select a combination of input features (traffic model,
topology, and queuing configuration).

To test the generalization capabilities to larger networks of our model, we
use a wide set of topologies of variable size (from 25 to 300 nodes). All
these topologies have been artificially generated using the Power-Law Out-
Degree Algorithm described in [14], where the ranges of the $\alpha$ and
$\beta$ parameters have been extrapolated from real-world topologies of the
Internet Topology Zoo repository [15].

In all the previous models, average traffic rates on src-dst flows are
carefully set to cover low to quite high congestion levels across different
samples, where the most congested samples have $\approx$$3\%$ of packet
losses.

###  IV-C Training

We implement the model using TensorFlow. To train the model, we use a custom
data augmentation approach that, given a link capacity ($Cap_{link}$), covers
a broad combination of $S_{f}$ and $C_{ref}$ values, in order to eventually
make the model generalize over samples with larger link capacities.
Particularly, given a link capacity, in some samples, we use low values of
$S_{f}$ with higher values of $C_{ref}$, while in other samples we make it in
the opposite way. As an illustrative example, if the model is trained over
samples with 1Gbps links, we can represent these capacities in different
samples as $Cap_{link}$=$10$*$100Mbps$$=$$1Gbps$, or
$Cap_{link}$=$1$*$1Gbps$$=$$1Gbps$. Thus, after training the model should be
able to make accurate inferences on samples that combine the maximum $S_{f}$
and $C_{ref}$ values seen during training – i.e.,
$Cap_{link}$=$10$*$1Gbps$$=$$10Gbps$. In practice, this means that the model
can be trained with samples with a maximum link capacity of 1 Gbps, and then
scale effectively to samples with link capacities up to 10 Gbps. Note that
these numbers are just illustrative, while this data augmentation method is
sufficiently general to produce in the training dataset wider ranges of
$S_{f}$ and $C_{ref}$ given a maximum link capacity. Thus, it can be
potentially exploited to represent combinations leading to arbitrarily larger
capacities.

Figure 4: Training and evaluation losses over time.

After making some grid search experiments, we set a size of 32 elements for
all the hidden state vectors ($h_{f}$, $h_{q}$, $h_{l}$), and $T$=$8$ message-
passing iterations. We implement $FRNN$, $LRNN$, and $U_{q}$ in Algorithm 1 as
Gated Recurrent Units (GRU) [16], and functions $R_{f}$ and $R_{q}$ as 2-layer
fully-connected neural networks with ReLU activation functions. Here, it is
important to note that the whole neural network architecture of the model
(Algorithm 1) constitutes a fully differentiable function, so it is possible
to train the model end-to-end. Hence, all the different functions that shape
its internal architecture are jointly optimized during training based on the
model’s inputs (i.e., the network state samples) and outputs (i.e., the
performance metrics).

We train the model for 200 epochs – with 4,000 samples per epoch – and set the
Mean Squared Error (MSE) as loss function, using an Adam optimizer with an
initial learning rate of 0.001. Figure 4 shows the evolution of the loss
during training on delay estimates (for the training and validation samples),
which shows stable learning along the whole training process.

Figure 5: Absolute relative error vs. topology size.

##  V Evaluation

As we have previously discussed, ML-based network models must generalize to
unseen and _larger_ networks to become a practical solution. In this vein, the
proposed GNN model was carefully designed to address this challenge. Also,
note that in the datasets provided for the Graph Neural Networking challenge
2021 all forwarding devices implement a FIFO queue scheduling policy, while
the proposed GNN-based model also supports different queueing policies with a
variable number of queues and priorities per port.

###  V-A Generalization to larger networks

In this section, we evaluate the model in a wide range of networks
considerably larger than the ones seen during training. Specifically, the
model was trained with topologies between 25 and 50 nodes and tested with
topologies from 50 to 300 nodes.

Figure 5 shows how the model generalizes to larger topologies not seen during
training. Particularly, the boxplots show the distribution of the absolute
relative error with respect to the topology size. As expected, the model
obtains better accuracy in topologies that are closer to the ones seen during
the training phase (50 to 99 nodes), achieving an average error of $4.5\%$
(green line). As the topology size increases, the average error stabilizes to
$\approx$$10\%$.

Generalization is an open challenge in the field of GNN. As we have previously
discussed in Sec. IV-A, we have addressed this by using domain-specific
knowledge and data augmentation. Particularly, we infer delay/jitter from
queues’ occupation and apply our scale-independent method to generalize to
larger topologies.

###  V-B Inference Speed

Finally, in this section, we evaluate the inference speed of the proposed GNN-
based model. Fast models are especially appealing for network control and
management, as they can be deployed in real-time scenarios. For this, we have
measured the execution times [Intel(R) Xeon(R) Gold 5220 CPU @ 2.20GHz] in the
experiments of the previous section. The results (Table I) show that the model
operate in the order of milliseconds. In particular, it goes from a few
milliseconds for small topologies to a few hundred for the larger ones.

| Topology Size
---|---
| [10,30] | [31,50] | [51,70]
Exec. Time | 48ms | 58ms | 100ms

TABLE I: Execution time vs. topology size

##  VI Related Work

The use of Deep Learning (DL) for network modeling has recently attracted the
interest of the networking community. This idea was initially suggested by
_Wang, et al._ [17]. The authors survey different techniques and discuss data-
driven models that can learn real networks. Initial attempts to instantiate
this idea use fully-connected neural networks (e.g. [18, 19]). Such early
attempts do not generalize to networks not seen in training, are not tested
with realistic traffic models, and do not model queues. More recent works
propose more elaborated neural network models, like Variational Auto-encoders
[20] or Convolutional NN [21]. However, they have similar limitations.

Finally, some early pioneering works use GNN in the field of computer networks
[22, 23]. However, these models consider a more simplified model of the
network, which do not account for different queuing policies, and the critical
property of generalizing to larger networks.

##  VII Discussion and Concluding remarks

In this paper, we have presented a novel GNN-based model that addresses a main
limitation of existing ML-based models applied to networks: generalizing
accurately to considerably larger networks ($\approx$10x) unseen during
training. Particularly, this solution would be suitable for the problem
proposed in the Graph Neural Networking challenge 2021 [2].

## References

  * [1] J. Suárez-Varela, M. Ferriol-Galmés, A. López, P. Almasan, G. Bernárdez _et al._ , “Must-read papers on gnn for communication networks,” https://github.com/BNN-UPC/GNNPapersCommNets, 2021.
  * [2] J. Suárez-Varela _et al._ , “The graph neural networking challenge: a worldwide competition for education in ai/ml for networks,” _ACM SIGCOMM Computer Communication Review_ , vol. 51, no. 3, pp. 9–16, 2021.
  * [3] F. Scarselli, M. Gori _et al._ , “The graph neural network model,” _IEEE Transactions on Neural Networks_ , vol. 20, no. 1, pp. 61–80, 2008\.
  * [4] P. W. Battaglia _et al._ , “Relational inductive biases, deep learning, and graph networks,” _arXiv preprint arXiv:1806.01261_ , 2018.
  * [5] P. Battaglia _et al._ , “Interaction networks for learning about objects, relations and physics,” in _Advances in neural information processing systems_ , 2016, pp. 4502–4510.
  * [6] D. Raposo _et al._ , “Discovering objects and their relations from entangled scene representations,” _arXiv preprint arXiv:1702.05068_ , 2017\.
  * [7] J. Gilmer _et al._ , “Neural message passing for quantum chemistry,” _arXiv preprint arXiv:1704.01212_ , 2017.
  * [8] M. Shreedhar and G. Varghese, “Efficient fair queuing using deficit round-robin,” _IEEE/ACM Transactions on networking_ , vol. 4, no. 3, pp. 375–385, 1996.
  * [9] J. Zhou _et al._ , “Graph neural networks: A review of methods and applications,” _arXiv preprint arXiv:1812.08434_ , 2018.
  * [10] L. Ruiz, L. Chamon, and A. Ribeiro, “Graph neural networks and the transferability of graph neural networks,” _Advances in Neural Information Processing Systems_ , vol. 33, 2020.
  * [11] L. Engstrom, B. Tran, D. Tsipras, L. Schmidt, and A. Madry, “Exploring the landscape of spatial robustness,” in _International Conference on Machine Learning_ , 2019, pp. 1802–1811.
  * [12] J. Su, D. V. Vargas, and K. Sakurai, “One pixel attack for fooling deep neural networks,” _IEEE Transactions on Evolutionary Computation_ , vol. 23, no. 5, pp. 828–841, 2019.
  * [13] A. Varga, “Discrete event simulation system,” in _European Simulation Multiconference (ESM)_ , 2001, pp. 1–7.
  * [14] C. R. Palmer and J. G. Steffan, “Generating network topologies that obey power laws,” in _IEEE Global Telecommunications Conference (GLOBECOM)_ , vol. 1, 2000, pp. 434–438.
  * [15] S. Knight, H. Nguyen _et al._ , “The internet topology zoo,” _IEEE Journal on Selected Areas in Communications_ , vol. 29, no. 9, pp. 1765 –1775, 2011.
  * [16] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” _arXiv preprint arXiv:1412.3555_ , 2014.
  * [17] M. Wang _et al._ , “Machine learning for networking: Workflow, advances and opportunities,” _IEEE Network_ , vol. 32, no. 2, pp. 92–99, 2017.
  * [18] A. Valadarsky _et al._ , “Learning to route,” in _ACM workshop on hot topics in networks_ , 2017, pp. 185–191.
  * [19] A. Mestres _et al._ , “Understanding the modeling of computer network delays using neural networks,” in _ACM SIGCOMM BigDAMA workshop_ , 2018, pp. 46–52.
  * [20] S. Xiao _et al._ , “Deep-q: Traffic-driven qos inference using deep generative network,” in _ACM SIGCOMM Workshop on Network Meets AI & ML_, 2018, pp. 67–73.
  * [21] X. Chen _et al._ , “Deep-rmsa: A deep-reinforcement-learning routing, modulation and spectrum assignment agent for elastic optical networks,” in _2018 Optical Fiber Communications Conference and Exposition (OFC)_ , 2018, pp. 1–3.
  * [22] F. Geyer and S. Bondorf, “Deeptma: Predicting effective contention models for network calculus using graph neural networks,” in _IEEE INFOCOM_ , 2019, pp. 1009–1017.
  * [23] K. Rusek, J. Suárez-Varela _et al._ , “Unveiling the potential of graph neural networks for network modeling and optimization in sdn,” in _ACM Symposium on SDN Research_ , 2019, pp. 140–151.