Bayesian Evidence Synthesis for Modeling SARS-CoV-2 Transmission

Anastasios Apsemidisa∗ and Nikolaos Demirisb

a∗Department of Statistics, Athens University of Economics and Business,
Athens, Greece<EMAIL_ADDRESS>
bDepartment of Statistics, Athens University of Economics and Business,
Athens, Greece<EMAIL_ADDRESS>

###### Abstract

The acute phase of the Covid-19 pandemic has made apparent the need for
decision support based upon accurate epidemic modeling. This process is
substantially hampered by under-reporting of cases and related data
incompleteness issues. In this article, a discrete-time stochastic epidemic
modeling framework is developed with discernible features targeted to publicly
available data. The models allow for estimating the total number of infections
while accounting for the endemic phase of the pandemic. We assess the
prediction of the infection rate utilizing mobility information, notably the
principal components of the mobility data. We elaborate upon vector analysis
of the epidemic dynamics, thus enriching the traditional tools used for
decision making. In particular, we show how certain 2-dimensional plots on the
phase plane may yield intuitive information regarding the speed and the type
of transmission dynamics.

Keywords: Covid-19; Epidemic model; Prediction; Variable selection; Principal
components analysis; Phase plane.

##  1 Introduction

The Covid-19 pandemic was initiated at Wuhan, China in late 2019 and is caused
by the spread of the SARS-CoV-2 virus. The exact burden of the pandemic
remains unknown and disease severity is highly variable with symptoms ranging
from asymptomatic to low fever to more serious such as chronic or death. The
high transmissibility of the disease (see Flaxman et al. 2020 for early
estimates) led to preventive measures being adopted, initially in the form of
non-pharmaceutical interventions. The need for monitoring the epidemic led to
the development of a vast range of models that could, in principle, assist the
decision making process on the effect of NPIs against the virus.

Traditional epidemic models are based on the recorded number of cases for
inference or prediction but in the case of Covid-19 it has become apparent
that only a proportion of the total cases is observed. This phenomenon occurs
due to the large number of asymptomatic cases as well as infected individuals
not being tested for a number of reasons. Thus, statisticians should build
upon these partly observed data to estimate quantities that, by definition,
depend on the total number of cases.

The standard deterministic Susceptible-Infectious-Removed (SIR) model is often
the basis of any analysis of infectious disease outbreak data and is typically
written as

$\begin{split}\frac{dS}{dt}&=-\lambda\cdot S\cdot I/N\\\
\frac{dI}{dt}&=\lambda\cdot S\cdot I/N-I/\tau\\\
\frac{dR}{dt}&=I/\tau\end{split}$ (1)

where $\lambda$ is the infection rate, $\tau$ is the infectious period, $S$,
$I$ and $R$ are the number of susceptible, infectious and removed individuals
respectively and $N$ is the size of the population. In this work we use a
suitably tailored stochastic discrete-time model.

The main contributions of this article are the following. A new discrete-time
stochastic epidemic framework is presented and parametrized based on publicly
available data. The model may be readily adapted to other infectious diseases
by amending the probability of death (defined below). The unobserved number of
infectious and susceptible individuals are estimated and independently
validated against external data. The marginal likelihood is analysed in detail
while a distinct type of variable selection procedure is used for predicting
the infection rate using mobility information. An alternative form of
Principal Components Analysis (PCA) is presented during this process, linking
the un-supervised nature of PCA to the supervised prediction task at hand. The
dynamics of the proposed system are examined in the phase-plane form and
suggested as a tool for monitoring the effectiveness of interventions as well
as a simple visual goodness of fit test for the model assumptions.

The remainder of the article is organised as follows. In Section 2 the models
considered in the article are presented along with methods for predicting the
infection rate and the total cases and an extension allowing for re-infection.
In Section 3, we present the results of training the models on data from the
United Kingdom (UK), Greece, and the United States of America (USA), the model
determination procedure and prediction performance. In Section 4 we present
the phase-plane analysis and its applications and we conclude with a
discussion in Section 5.

##  2 Modeling framework

In this section we describe the proposed methodology for epidemic modeling at
the country level, a suitably tailored stochastic SEIR in discrete time. The
model comprises of four states, namely: Susceptible, Exposed, Infectious and
Removed, abbreviated as $S$-state, $E$-state, $I$-state and $R$-state
respectively. Note that while we focus upon the application to the Covid-19
data the model is applicable to other communicable diseases.

###  2.1 Stochastic discrete-time transmission model

Let $\mathbf{D}$ be a random vector of daily deaths in a given country and $d$
its recorded realization. Let also $\mathcal{F}_{t}$ be the history of the
process $\\{D_{t},\theta_{t}\,|\,t\in\mathbb{N}\smallsetminus\\{1\\}\\}$ up to
time $t$, where $\bm{\theta}$ is the vector of the conditional means of the
process, i.e. $\theta_{t}=\mathbb{E}[D_{t}\,|\,\mathcal{F}_{t-1}]$. We model
the mean number of daily deaths using a Negative Binomial distribution with
parametrization for a random variable $X$ and realizations $x$

$NB(x;\mu,\psi)=\frac{\Gamma(x+\psi)}{x!\,\Gamma(\psi)}\left(\frac{\mu}{\mu+\psi}\right)^{x}\left(\frac{\psi}{\mu+\psi}\right)^{\psi}$

where $\mathbb{E}[X]=\mu$ and $Var[X]=\mu+\displaystyle\frac{\mu^{2}}{\psi}$
with $(\mu,\psi)\in\mathbb{R}\times\mathbb{R}^{+}$. Thus, the likelihood reads

$D_{t}\,|\,\mathcal{F}_{t-1}\sim NB(\theta_{t},\psi)$ (2)

We wish to make inferences about the disease transmissibility based upon the
total number of cases, as opposed to the observed ones. Thus, we assume that
the mean number of deaths is a proportion of the total cases that were
generated during the previous days. We initiate the model assuming an equal
number of cases, $C$, for the first $\tau+h$ days, i.e. $C_{t}=C_{1}$, for
$t=1,...,\tau+h$. Then, for a sample of $n$ days the model may be written as

$\displaystyle\theta_{t}$
$\displaystyle=p_{t}\cdot\sum_{k=1}^{t-1}\pi_{t-k}\cdot
C_{k}\,,\quad\mathrm{for}\quad t=2,...,n$ (3) $\displaystyle C_{t}$
$\displaystyle=\lambda_{t-1-h}\frac{S_{t-1-h}I_{t-1-h}}{N}\,,\quad\mathrm{for}\quad
t=\tau+h+1,...,n-1$ (4)

where $p_{t}\in(0,1)$ is the Infection Fatality Ratio (IFR) at time $t$, i.e.
the probability of death for a given case that occurred at time $s<t$,
$C_{t}\in\mathbb{R}^{+}$ are the total cases at time $t$,
$\lambda_{t}\in\mathbb{R}^{+}$ denotes the infection rate at time $t$,
$S_{t}\in\mathbb{R}^{+}$ and $I_{t}\in\mathbb{R}^{+}$ are the numbers of
susceptible and infectious individuals respectively at time $t$, while
$N\in\mathbb{N}$ is the (fixed) population size of the country under study.
The length of the Exposed period, during which an individual is infected but
not infectious, is represented by $h\in\mathbb{N}_{0}$, while
$\tau\in\mathbb{N}$ is the infectious period. Both of those time periods are
assumed to be known from previous studies.

For the infection fatality ratio, we assume a piecewise constant form with
$B-1$ change-points at times $l_{b}$, i.e. $p_{t}=p_{(b+1)}\cdot
I\big{(}t\in[l_{b},l_{b+1}-1]\big{)}$, for $b=0,...,B-1$ with $l_{0}=2$,
$l_{B}=n$ and $p_{n}=p_{(B-1)}$ (note that $l_{0}$ and $l_{B}$ are not actual
change-points, but they are used to define the intervals of constant IFR’s).
For the $\pi_{t-k}$ term that posits different weights to each past day, we
use the sum of two Gamma distributions $\pi(t)$ as an estimate of the time
from infection until death and discretize it as

$\pi_{s}=\int\limits_{s-0.5}^{s+0.5}\pi(t)dt\quad,\,\mathrm{for}\quad
s=2,...,n-1$

and $\pi_{1}=\int\limits_{0}^{1.5}\pi(t)dt$. We use these weights $\pi_{s}$ in
the model as known proportions. The infection rate is written as a piecewise
constant function of time with $J-1$ change-points at times $u_{j}$, i.e.
$\lambda_{t}=\lambda_{(j+1)}\cdot I\big{(}t\in[u_{j},u_{j+1}-1]\big{)}$, for
$j=0,...,J-1$ with $u_{0}=1$ and $u_{J}=n-h-1$ (like with the IFR case,
$u_{0}$ and $u_{J}$ are not change-points). We update at every time step the
number of susceptible individuals $S_{t}$ and the active set of infectious
individuals $I_{t}$ as

$\displaystyle S_{t}$ $\displaystyle=S_{t-1}-C_{t}$ (5) $\displaystyle I_{t}$
$\displaystyle=\sum_{k=0}^{\tau-1}C_{t-k}$ (6)

for $t=\tau,...,n-h-2$. An important quantity in epidemiology of infectious
diseases is the basic reproduction number $R_{0}$, typically interpreted as
the number of secondary infections generated by an infectious individual in a
large susceptible population. If the reproduction number $R_{t}$ at day $t$ is
greater than 1 (less than) the epidemic is increasing (decreasing). This
quantity is estimated using $R_{t}=\lambda_{t}\cdot\tau\cdot S_{t}/N$.

The description of the basic model is completed by the $R$-state, estimating
the number of individuals that died or recovered by summing those leaving the
Infectious state, i.e.

$R_{t}^{(s)}=\sum_{i=1}^{t-\tau}C_{i}$ (7)

for $t=\tau,...,n-h-2$. The superscript “$(s)$” denotes that this is the
Removed state at time $t$, as opposed to the reproduction number $R_{t}$.

Equations (2) through (7) define the dynamic stochastic discrete-time SEIR
model. Note also that setting $h=0$, the model collapses to the simple SIR
type. The model represents a natural stochastic discrete time analog of the
deterministic SIR system given by the ordinary differential equations (ODE) at
(1).

####  2.1.1 Incorporating vaccination and demography

A turning point of the pandemic was the introduction of the vaccine, which
played a vital role in mitigating its influence. The vaccinations are added in
our model by assuming that after the first dose each individual remained
susceptible for two weeks and then they became immune with probability
$a_{1}$. After three more weeks (corresponding to roughly the time the second
dose was done in most countries), we assume that this probability of immunity
is being raised to $a_{1}+a_{2}$. Thus, we incorporate the vaccine data by
directly moving individuals to the $R$-state with a fixed probability. The
updates of the susceptible equation (5) are now performed by

$S_{t}=S_{t-1}-C_{t}-V_{t}$ (8)

where

$V_{t}=a_{1}\cdot\rho_{t-14}\cdot
I\big{(}t\in[15,n-h-2]\big{)}+a_{2}\cdot\rho_{t-35}\cdot
I\big{(}t\in[36,n-h-2]\big{)}$

and $\rho_{t}\in\mathbb{N}_{0}$ is the number of vaccinations at time $t$. The
removed individuals at time $t$ are now
$R_{t}^{(s)}=\displaystyle\sum_{i=1}^{t-\tau}C_{i}+V_{t}$, while the
infectious remain the same as in (6). We add further realism to the model by
the inclusion of demography. This is not necessary when modeling acute
outbreaks of short duration but over a period of three years one may wish to
consider including births and deaths into the population dynamics. We assume
that the number of births $A\in\mathbb{R}^{+}$ equal the one of deaths (due to
reasons other than coronavirus) and update equation (8) as follows:

$S_{t}=S_{t-1}-C_{t}-V_{t}+A\cdot(1-S_{t-1}/N)$ (9)

The $A$ term accounts for new births while the term $A\cdot S_{t-1}/N$
accounts for deaths. Furthermore, deaths from natural causes can occur inside
the active set, so $I_{t}$ is updated by

$I_{t}=\sum_{k=0}^{\tau-1}C_{t-k}-A\cdot I_{t-1}/N$ (10)

while the removed population is given by

$R_{t}^{(s)}=\sum_{i=1}^{t-\tau}C_{i}+V_{t}-A\cdot R_{t-1}^{(s)}/N$ (11)

Note that newborn individuals are assumed to be susceptible and so we do not
add births to $I_{t}$ or $R_{t}^{(s)}$. Thus, equations (2), (3), (4), (9),
(10) and (11) synthesize the SEIR model with vaccination and demography.

The $A$ term is assumed to be a known constant informed by the data while
$V_{t}$ is simply a transformation of data. The fixed parameter $A$ can be
estimated using the population distribution across the published age groups as
follows. Suppose that there exist $K$ age groups of length $g_{i}$ years, for
$i=1,...,K$, and in each group belong $A_{i}$ people. Then, we assume that

$A=\frac{A_{1}}{365\cdot g_{1}}$ (12)

are the newly born people set equal to the number of daily demographic deaths.

When adding vaccinations and demography to the model, $\mathcal{F}_{t}$ is the
$\sigma$-algebra generated by
$\\{D_{t},\theta_{t},\mathbf{X}_{t}\,|\,t\in\mathbb{N}\smallsetminus\\{1\\}\\}$,
where $\mathbf{X}_{t}$ is the “covariate” information of the weights $\pi_{t}$
and vaccinations $V_{t}$. Note that equation (2) refers to the distribution of
deaths conditional upon the parameters $\psi$, $\lambda_{t}$, $p_{t}$, $C_{1}$
and $\theta_{1}$, so formally it should read

$D_{t}\,|\,\psi,\lambda_{t},p_{t},C_{1},\theta_{1},\mathcal{F}_{t-1}\sim
NB(\theta_{t},\psi)$ (13)

####  2.1.2 Prediction of the infection rate

The SEIR model we have introduced assumes a piece-wise constant infection rate
$\lambda$, which in principle can be affected by certain covariates. A crucial
such factor possibly determining $\lambda_{t}$ is the mobility of the
population, since the rate of infection is typically thought of as the product
of the probability of infection given there is a contact times the number of
contacts made by an individual. The mobility data provided by Google  and
Apple  during the course of the epidemic are daily percent changes of mobility
in certain places compared to a baseline mobility level. Google records six
variables, while Apple records two, all of which are correlated to each other.
Therefore, some form of de-correlation is required. To this end we used a
Principal Components Analysis (PCA) since it also facilitates dimension
reduction.

We initially embedded this process within the wider model and used the first
Principal Component (PC) of the these data as the mobility variable $m_{t}$ at
day $t$. We modelled the infection rate as a piecewise constant Autoregressive
Moving Average (ARMA) process, i.e.

$\mathrm{log}(\lambda_{t})=\bigg{\\{}\beta_{0j}+\beta_{1j}m_{t}+\sum_{m=1}^{M}\phi_{mj}\mathrm{log}(\lambda_{t-m})+\sum_{k=1}^{K}\delta_{kj}d_{t-k}\bigg{\\}}\cdot
I\big{(}t\in[u_{j},u_{j+1}-1]\big{)}$

for $j=0,...,J-1$. Overall, it appears that the information in the data was
insufficient to estimate the parameters that are several levels away from the
data. In particular, this phenomenon has manifested itself by convergence
problems of the algorithm.

We investigated the post-processing option of first fitting the SEIR model
with a piecewise constant infection rate and then treating the produced mean
$\lambda_{t}$ estimates as a response variable in a regression scenario, i.e.
ignoring the uncertainty of the $\lambda_{t}$ estimate. In particular, we
assumed that $\mathbb{E}[\mathrm{log}(\hat{\lambda}_{t})|m_{t}]=g(m_{t})$,
where $g(\cdot)$ is a suitable function. We investigated three types of
predictive models, namely Linear Regression (LR), Generalized Additive Models
(GAM) and Extreme Gradient Boosted Regression Trees (XGB) and compared them in
terms of Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE)
estimated by $r$-repeated $k$-fold Cross Validation (CV). We consider three
scenarios of constructing the covariates of mobility: working with the first
PC of mobility, the first two PCs, or using a specific variable chosen as
described next. For each of these scenarios we test whether to smooth the
variable(s) $m_{t}$ or not via the discretized serial interval distribution
$f_{t}$ of Covid-19. In the case of smoothing we work with
$m_{t}^{*}=\displaystyle\sum_{i=1}^{t-1}f_{t-i}m_{i}$, for $t=2,...,n-2-h$ and
$m_{1}^{*}=m_{1}$. The discretization is performed the same way as with the
$\pi_{s}$ case. Finally, when working with the PCs we further considered two
scenarios: including in the PCA procedure all of the mobility data or
excluding the one referring to the places of residence. The intuition is that
when this variable is excluded $m_{t}$ has the natural mobility meaning, in
the sense that increasing $m_{t}$ leads to increased actual mobility.
Otherwise, the linear combination of the mobility variables would include a
variable that does not correspond to mobility, thus perplexing the
interpretation.

Variable selection is performed in a way that allows the mobility variable
used to contain information for both the omitted covariates and the response.
Adopting matrix notation we denote $m_{ik}$ the value of the variable $k$ at
sample $i$, i.e. $\mathbf{m}_{k}$ (for $k=1,...,8$) is one of the variables to
be selected. Let $\mathrm{log}(\lambda)=\mathbf{y}$ and
$\hat{\mathbf{y}}^{(k)}$ be the fitted values of the linear regression of
$\mathbf{y}$ versus $\mathbf{m}_{k}$, while $\hat{\mathbf{m}}^{(k)}$ be the
fitted values of the multiple linear regression of $\mathbf{m}_{k}$ versus all
the other mobility variables except $\mathbf{m}_{k}$. Then, we seek to solve
the following problem

$\operatorname*{arg\,max}_{k}\quad\Bigg{\\{}\,\frac{\displaystyle\sum_{i=1}^{n-2-h}(\hat{m}_{i}^{(k)}-\bar{m}_{k})^{2}}{\displaystyle\sum_{i=1}^{n-2-h}(m_{ik}-\bar{m}_{k})^{2}}+\frac{\displaystyle\sum_{i=1}^{n-2-h}(\hat{y}_{i}^{(k)}-\bar{y})^{2}}{\displaystyle\sum_{i=1}^{n-2-h}(y_{i}-\bar{y})^{2}}\,\Bigg{\\}}$
(14)

where the bar over a variable indicates the mean over the sample data.

Moreover, we propose three new measures of the induced distortion of
discarding $k$ out of the original $L$ variables measured on $n$ observations.
Let $\mathbf{X}$ be the $n\times L$ matrix of the centered and scaled data.
PCA applies a transformation $\mathbf{A}$ on $\mathbf{X}$ (the eigen-
decomposition of the covariance) resulting in a new data matrix
$\mathbf{\tilde{X}}$. In the case of no dimension reduction
$\mathbf{X}\cdot\mathbf{A}=\mathbf{\tilde{X}}$ is simply a rotation of the
original data, so applying the inverse results again in $\mathbf{X}$.
Otherwise $\mathbf{\tilde{X}}\cdot\mathbf{A}_{k}^{-1}$ represents a distorted
image of $\mathbf{X}$. Let $\mathbf{\tilde{X}}_{k}$ be the distorted dataset
with $k$ the number of retained basis vectors. Note that
$SS_{k,l}=(\mathbf{X}_{l}-\mathbf{\tilde{X}}_{k,l})^{2}\cdot\mathbf{1}^{\prime}$
is a similarity measure between the original variable $l$ of $\mathbf{X}$, and
the new variable $l$, $\mathbf{\tilde{X}}_{k,l}$ quantifying the information
lost from variable $l$ when $k$ out of $L$ variables are retained
($\mathbf{1}$ is a vector of $n$ 1’s). Thus, the decision of keeping $k$ or
$k-1$ components can be assisted by the information lost from a specific
variable calculating $(SS_{k,l}-SS_{k-1,l})/SS_{k-1,l}$. Generalizing, we can
find the percent change on the distortion over all variables by
$\sum_{l=1}^{L}(SS_{k,l}-SS_{k-1,l})/\sum_{l=1}^{L}SS_{k-1,l}$, where
$SS_{k,l}$ is the distortion of variable $l$. This kind of information is
already given by the proportion of variance explained, although in a different
form.

However, the distortion of each variable also helps in scenarios where the
number of components should be decided based on adequate participation of a
specific variable in the resulted PC’s, for instance when another analysis is
to be performed based on the reduced dataset. It moves PCA from a purely
unsupervised technique towards a more supervised scheme without changing the
algorithm itself like for instance the partial least squares idea. Lastly, in
$SS_{k,l}$ we sum the squared differences of the old and new data for variable
$l$, but any other similarity measure can be used like in a clustering
scenario (maximum distance, mean distance, minimum distance etc). For
instance, instead of the sum of squares, we can measure the distortion of each
variable via $D_{l}=\max_{i}\,|(\mathbf{X}_{l}-\mathbf{\tilde{X}}_{k,l})|$ for
$i=1,...,n$.

Another natural way to quantify the distortion of each variable is by the
magnitude of the scatter induced by the transformation $\mathbf{A}_{k}^{-1}$.
Since $\mathbf{A}_{L}^{-1}$ returns the original data, plotting
$\mathbf{X}_{l}$ against $\mathbf{\tilde{X}}_{k,l}$ will result in points on
the $y=x$ line only when $k=L$. Else, there will be a deviation which can be
measured either by comparing the sum of squares or by the angle between the
fitted line of their regression and $y=x$. Thus, two measures of distortion
for a specific variable $l$ are

$D_{1}=SS_{k,l}^{-1}\cdot(\mathbf{\tilde{X}}_{k,l}-\mathbf{\hat{\tilde{X}}}_{k,l})^{2}\cdot\mathbf{1}^{\prime}$

and

$D_{2}=|\mathrm{Arctan}(1)-\mathrm{Arctan}(\hat{b})|$

where $\mathbf{1}$ is a vector of 1’s with length $n$ and $\hat{b}$ is the
slope estimated from the aforementioned regression. Which of the above should
be used is not only a matter of personal preference, but also of their
sensitivity to the variables under study, in the sense that one measure may
indicate a larger distortion relatively to the others. Their relative
strictness and variability is currently under study. Results of the PCA
procedure are given in the Appendix.

####  2.1.3 Prediction of the Total Cases

The number of tests performed each day affects the number of recorded Covid-19
cases. We assess this effect as follows: we aim to predict the total number of
cases $C_{t}$ using the number of observed cases $c_{t}$, deaths $d_{t}$,
tests $T_{t}$, intensive care unit (ICU) entries $H_{t}$ and the first two
PC’s of the 7 mobility variables smoothed by the serial interval distribution
(the places of residence variable is not included). The three competing
predictors are again LR, GAM and XGB.

###  2.2 Waning immunity and the SEIRS model

The SEIR model presented thus far does not allow a transition from the
$R$-state to the $S$-state. While a reasonable approximation at the early
phase of the pandemic, re-infections cannot be ignored for longer time-
horizons. Thus, we extend the SEIR model with vaccination and demography to
the SEIRS setting where recovered individuals move to susceptibility after
$t^{*}$ days using $r_{t}=(1-p_{t})\cdot\sum_{k=1}^{t-1}\pi_{t-k}^{*}\cdot
C_{k}$, where $\pi_{t-k}^{*}$ is the discretized Gamma distribution of time
from infection until recovery estimated in Paul and Lorin (2021). Adding the
lagged $r_{t}$ to $S_{t}$ in equation 9 via

$S_{t}=S_{t-1}-C_{t}-V_{t}+A\cdot(1-S_{t-1}/N)+r_{t-t^{*}}$ (15)

allows recovered individuals to lose their immunity $t^{*}$ days after
infection.

###  2.3 Bayesian inference

####  2.3.1 Computation

The SEIR model described above is fitted under the Bayesian paradigm. The
complexity of this non-linear hierarchical stochastic model implies that no
analytic likelihood calculations exist. We can use Hamiltonian Monte Carlo
(HMC) sampling for posterior inference (see Appendix S1), notably the NUTS
algorithm (see Hoffman et al. 2014), which is implemented in the Stan
probabilistic programming language.

####  2.3.2 Prior specification and elicitation

We assume a priori independent priors on the $B$ infection fatality ratios
$p_{(b)}$, the $J$ infection rates $\lambda_{(j)}$, the dispersion parameter
$\psi$ and the initial value of total cases $C_{1}$ and fitted deaths $d_{1}$.
The first $\tau+h$ values of $C_{t}$ are assumed equal and are given a
$Gamma(2,0.0625)$ prior. The initial value of the susceptible population is
then set to $S_{1}=N-C_{1}$. For each lambda between change-points $u_{j}$ and
$u_{j+1}-1$ we assign a $LogNormal(0,1)$ prior. The dispersion parameter
$\psi$ of the Negative Binomial is allocated a $Gamma(2,0.125)$ prior. While
these prior distributions are weakly informative, we also conducted
sensitivity analyses to assess their effect as typically done in Bayesian
robustness settings. We use a fixed $\tau=6$ infectious period (Cereda et al.
2020; Flaxman et al. 2020), and a fixed exposed period of $h=2$, since the
mean incubation period is approximately 5 days (Lauer et al. 2020) and
infectiousness starts approximately 2 days before the symptom onset (He et al.
2020).

The prior on the IFR parameters $p_{(b)}$ requires particular caution. For
each interval of constant IFR, the prior was set to a strongly informative
Gaussian distribution with standard deviation of $10^{-4}$ and mean computed
as follows. First, we scaled the IFR published by the Centers for Disease
Control and Prevention (CDC) for the 4 age groups 0-17, 18-39, 40-64 and 65+,
$p^{(0)}_{k}$ ($k=1,...,4$), according to the age distribution of the recorded
cases in the country under study. Then, we locate the change-points by
inspecting the observed country-specific IFR and use the mean IFR inside each
time interval. Thus, if we denote by $c_{t,k}$ the cases at time $t$ for age
group $k$ in a given country, the IFR is computed as
$p_{t}^{*}=\sum_{k=1}^{4}p^{(0)}_{k}\displaystyle\frac{c_{t,k}}{\sum_{i=1}^{4}c_{t,i}}$
and we set the mean of the Gaussian prior to
$\mathbb{E}[p_{(b)}]=\displaystyle\frac{1}{l_{b+1}-l_{b}}\sum_{t=l_{b}}^{l_{b+1}-1}p_{t}^{*}$.

##  3 Results

In this section, we analyse the pandemic data from the UK, Greece and the USA.
Modeling the entire USA as a single entity may not be entirely appropriate due
to the inherent heterogeneity stemming from both the large geographic area and
population size so these results should be interpreted with caution.

###  3.1 Sources of evidence

The publicly available data we leverage include the recorded cases and deaths,
the number of vaccinated individuals, the age distribution of cases, the
demographic births and deaths, the country population, the number of tests
performed and quantities estimated at the beginning of the pandemic, which are
the time from infection until death, the infectious period, the exposed period
and the serial interval. Below we outline the values we use for these
quantities (also see Appendix S2).

For the UK the cases and deaths data as well as the vaccination and age
distribution of cases are provided in the Public Health England  website. For
the parameter $A$ we use demographic evidence  and calculate it using
$g_{1}=5$ and $A_{1}=\displaystyle\frac{6.2\cdot N_{UK}}{100}$. We use
$N_{UK}=67886011$ as the 2020 UK population size .

The main component of the likelihood function is linked to the daily number of
deaths. For Greece these are obtained from the COVID-19 Data Repository of the
Center for Systems Science and Engineering at Johns Hopkins University (Dong
et al. 2020). The vaccine doses received are taken from Our World in Data
(Mathieu et al. 2021), while the age distribution of cases and the daily tests
performed were taken by the Github page of Sandbird . The demographic births
and deaths are obtained from the Hellenic Statistical Authority (HSA ) which
provides the population age groups and so we use equation (12) to calculate
$A$ with $g_{1}=10$ and $A_{1}=1049839$. The population of Greece in 2020 is
set to $N_{GR}=10816286$ following the HSA estimate. For the USA deaths data
are taken from the Johns Hopkins University while the cases and the
vaccinations data are obtained from Our World in Data. The population is set
to be $N_{USA}=331230000$ and we calculate $A$ using $g_{1}=15$ and
$A_{1}=\displaystyle\frac{18.51\cdot N_{USA}}{100}$. For these we use the
statista.com page, specifically the published total population  and age
distribution  in 2020.

For the infection-to-death distribution we use the sum of the infection-to-
onset of symptoms and onset-to-death distributions, which are Gamma with
shapes 1.35 and 4.94 and rates 0.27 and 0.26 respectively, as in (Flaxman et
al. 2020). The same source was used for the serial interval distribution as
well. The $p_{k}^{(0)}$ IFR can be found in the CDC  page. The probabilities
of moving to the $R$-state are set to be $a_{1}=0.4$ and $a_{2}=0.1$.

###  3.2 Epidemic parameters and functions thereof

Figure 1: Posterior estimates of $R_{t}$ (left) and $\theta_{t}$ (right) for
Greece during the first two years of the epidemic. The median is depicted with
a solid black line, along with 50% and 95% credible intervals.

We focus on the acute phase of the pandemic and fit the model to data from
Greece covering the 26/2/2020 to 31/12/2021 period. We independently validate
our estimates based on a large UK seroprevalence study. The NUTS algorithm is
mostly used to estimate model parameters. We also used Simulated Annealing
(SA) in order to maximise the un-normalized posterior. Extensive searches by
SA required less than 10 minutes when sampling took more than two days.
However, the SA results were very sensitive to initial values while HMC
appeared stable and robust. We therefore report HMC-based results due to
statistical but not computational efficiency.

Figure 2: Left: Cumulative total cases: posterior median (solid black line)
with 50% and 95% credible intervals. Right: Proportion of observed cases
smoothed by local regression.

In Figure 1 the estimated reproduction number, $R_{t}$, and mean deaths
$\theta_{t}$ are depicted. The reasonable fit to the death data is reassuring.
The piecewise constant $R_{t}$ is scaled by the proportion of susceptible
individuals as the theory suggests. The total number of cases at day $t$,
$C_{t}$, i.e. is the sum of the recorded and unrecorded ones are depicted in
the left plot of Figure 2 while the right panel contains the smoothed
estimated proportion of observed cases, $c_{t}/C_{t}$. We estimate that the
first million infections (10% of the total population) was reached on April
2021 but observed eight months later on December 2021. The probability of
recording a case initially was around 1/4 but then increased as tests and
self-tests became widely available, closer to 3/4 of the total cases.

Figure 3: Cumulative total cases for the UK. The time when the REACT-2 study
estimated the total cases is noted as a vertical line, while the reported 95%
confidence interval is shown as horizontal dashed lines.

The external validity of our model was independently as assessed by training
the model on the UK death data and superimposing on Figure 3 the estimated
number of infections against the REACT study (Ward et al. 2021) which
estimated an approximate 6% prevalence on July 2020. It is apparent that our
model estimate agrees well with this independent data source.

###  3.3 Sensitivity Analysis and Model Comparison

As part of the model determination process we tested the initial values, the
prior distributions, we added vaccination and demography to the model and
tested seven scenarios where we removing the exposed period, vaccination or
demography. In Table 1 we compare eight models: four SIR and four SEIR each
with zero to two of the extensions “vaccination” and “demography”. The
comparison is made using the information criteria AIC, BIC, DIC, DIC using
half the variance of the deviance (referred to as DIC2) and WAIC. The time
needed to fit each model is also noted. In Appendix S3, we also compare the
logarithm of the evidence $p(d_{t})$ for each model using Bridge sampling and
pairwise Bayes factors. Overall, the SEIR with vaccination is selected by AIC,
BIC and DIC, the SEIR with demography model by DIC2, while the SEIR with
vaccination and demography is preferred by WAIC.

All the training times are similar and approximately equal to 2.5 days. The
SEIR model with vaccination has the largest marginal likelihood, slightly
higher than that of the SEIR with vaccination and demography while the
interval estimates of the corresponding $\log p(d_{t})$ show substantial
overlap. Fitting the SEIRS model to the Greek data, allowing for re-infection
3 months after recovery ($3\cdot 4\cdot 7$ days) did not make any material
difference to the SEIR version except for the last month, when the
characteristics of the epidemic also change.

Table 1: Information criteria for the eight tested models and their training time measured in days. The results are rounded to 2 decimal digits. | AIC | BIC | DIC | DIC2 | WAIC | Time $(d)$
---|---|---|---|---|---|---
SIR | 4034.75 | 4283.06 | 3924.75 | 3974.89 | 3980.82 | 2.67
SIR.vacc | 4034.66 | 4282.97 | 3924.66 | 3976.76 | 3980.66 | 2.43
SIR.dem | 4034.93 | 4283.24 | 3924.93 | 3978.20 | 3981.05 | 2.57
SIR.vacc.dem | 4035.06 | 4283.37 | 3925.06 | 3977.96 | 3981.06 | 2.54
SEIR | 4034.91 | 4283.22 | 3924.91 | 3976.59 | 3979.68 | 2.39
SEIR.vacc | 4034.62 | 4282.93 | 3924.62 | 3974.38 | 3979.09 | 2.37
SEIR.dem | 4034.52 | 4282.83 | 3924.52 | 3974.88 | 3978.61 | 2.50
SEIR.vacc.dem | 4034.56 | 4282.87 | 3924.56 | 3975.72 | 3977.84 | 2.39

The distributional form of the likelihood of the death data was tested as
follows. The Negative Binomial is a scale mixture distribution of a Poisson
with Gamma mixing rate and we also examined the standard Poisson distribution
as well as scale mixtures with comparable Exponential and Log-Normal options.
Thus, this part of the model now reads $D_{t}\sim P(\theta_{t})$ with
$\theta_{t}\sim Exp(1/\mu_{t})$ or $\theta_{t}\sim LN(m_{t},s_{t}^{2})$, where
$m_{t}=\displaystyle\frac{\log(\mu_{t}^{2})}{\displaystyle\sqrt{(}\mu_{t}^{2}+\sigma^{2})}$
and $s_{t}^{2}=\log(1+\sigma^{2}/\mu_{t}^{2})$. Then,
$\mu_{t}=p_{t}\sum_{k=1}^{t-1}\pi_{t-k}C_{k}$ as we have with the Negative
Binomial case. The Poisson-Exponential model (a special case of the Negative
Binomial) fits the data well and reduces training time to a half. On the other
hand, the Poisson-LogNormal struggled to converge and was unreliable.

###  3.4 Prediction of $\lambda_{t}$ and $C_{t}$

We used PCA to reduce the 8 dependent mobility variables to 1 or 2 keeping
78.09% or 89.25% of the initial variance. The prediction of the infection rate
using mobility data is summarised in Table (2). Both RMSE and MAE suggest that
the best prediction is achieved by the XGB model using the 2 smoothed PC’s
where the places of residence are excluded. The overall prediction ability is
generally low. The mobility variable that describes best all the others and
simultaneously serves as a good predictor of the infection rate is the one
associated with transit stations, which is rather intuitive since public
transport relates well to contacts (see Table 2). The XGB model outperforms
its competitors in predicting total cases too (see Table 3).

Broadly, the first PC summarizes the actual mobility in the sense that it
takes positive values when mobility increases, while negative values
correspond to increased stay in places of residence. The second PC is used
only for prediction purposes and does not have a clear interpretation. Thus,
fitting a mixture of two Gaussian distributions using the Expectation-
Maximization algorithm to the scores of the first PC, we find that “mobility”
mainly stems from $0.31\cdot N(-2.86,1.2^{2})+0.69\cdot N(1.26,1.78^{2})$. In
Appendix S4, we visualize the two Gaussian components that prevailed at each
day during Covid-19 waves in Greece and we can see for instance that during
the two lock-downs in spring 2020 and winter of 2020 to 2021, mobility was
mostly generated by the first component of the mixture, which corresponds to
the prevalence of the places of residence.

Table 2: Prediction of the infection rate: 5-repeated 10-fold CV prediction error estimates for each model and scenario as captured by RMSE and MAE (in square brackets). |  |  | LR | GAM | XGB
---|---|---|---|---|---
1 PC | Ser Int | All in | 0.3267 [0.2620] | 0.3120 [0.2484] | 0.3043 [0.2410]
No resid | 0.3284 [0.2631] | 0.3140 [0.2498] | 0.2998 [0.2372]
No Ser Int | All in | 0.3336 [0.2723] | 0.3312 [0.2698] | 0.3197 [0.2624]
No resid | 0.3346 [0.2734] | 0.3297 [0.2685] | 0.3178 [0.2607]
2 PC | Ser Int | All in | 0.3153 [0.2500] | 0.2900 [0.2220] | 0.2800 [0.2142]
No resid | 0.3128 [0.2479] | 0.2927 [0.2260] | 0.2699 [0.2064]
No Ser Int | All in | 0.3333 [0.2712] | 0.3267 [0.2662] | 0.3162 [0.2595]
No resid | 0.3337 [0.2713] | 0.3247 [0.2646] | 0.3137 [0.2563]
Variable Selection |  | Ser Int | 0.3231 [0.2608] | 0.3072 [0.2365] | 0.3018 [0.2381]
| No Ser Int | 0.3299 [0.2699] | 0.3202 [0.2654] | 0.3128 [0.2597]
Table 3: Prediction of the total cases: 5-repeated 10-fold CV prediction error estimates for each model and scenario as captured by RMSE and MAE. | RMSE | MAE
---|---|---
LR | 1.0741 | 0.8309
GAM | 0.5295 | 0.3888
XGB | 0.5040 | 0.3669

###  3.5 Transition to endemicity

Figure 4: Daily estimates of the covariance (left) and correlation (right)
between $\lambda_{t}$ and $\tau S_{t}/N$ for USA.

Defining the end of the acute phase of the epidemic is a non-standard problem.
The definition of the reproduction number $R_{t}=\lambda_{t}\cdot\tau\cdot
S_{t}/N$ implies that at the early stage of the epidemic $R_{t}$ is
approximately proportional to the infection rate. The covariance of the
infection rate and the scaling factor $\tau S_{t}/N$ offers an insight to when
those two diverge (see Figure 4) whence the covariance becomes negative. An
additional version of the divergence between $\lambda_{t}$ and $R_{t}$ may be
found in Appendix S5.

##  4 Phase Plane Analysis

In this Section we examine the model behaviour through the lens of dynamical
non-linear systems. In particular, the non-pharmaceutical interventions are
expressed in a novel way and a qualitative analysis of the epidemic is offered
through the visual inspection and interpretation of certain plots. We
introduce the main concepts in continuous time, before adapting to the
discrete time used in this paper.

###  4.1 Epidemics on the Continuous Time Domain

The system of coupled equations of $S$ and $I$ may be inspected focusing on
the bivariate space $S\times I$ and the dynamic behaviour on this $SI$-plane,
denoted as the set $V\subseteq\mathbb{R}_{+}^{2}$. The system (1) can be
examined as a vector field $\mathbf{F}$ on $V$. We define the natural epidemic
flow to be a curve $\bm{\sigma}(t)=\big{(}S(t),I(t)\big{)}$ that satisfies

$\mathbf{F\big{(}\bm{\sigma}(t)\big{)}}=\frac{dS}{dt}\cdot\bm{\hat{\textbf{\i}}}+\frac{dI}{dt}\cdot\bm{\hat{\textbf{\j}}}$

Given an initial condition $\bm{\sigma}(0)$ the epidemic flow defines a
trajectory on $V$, which we call the natural course of the epidemic. The image
of the realized transformation of time $t$ to the $\mathbb{R}_{+}^{2}$ plane
$\bm{\sigma}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}^{2}$ (that is the
realized trajectory) is a useful plane curve since it defines a path whose $x$
and $y$ coordinates represent the number of susceptible and infected
individuals respectively at each time $t$. Thus, monitoring the path a
particle follows in the epidemic flow, starting from the initial position
$(S(0),I(0))$ can be informative.

The actual behaviour of the epidemic model of the present paper differ since
the parameters change and intervention measures interrupt the theoretical
dynamics. The definition of the actual course of the epidemic relates to a
curve $\bm{\gamma}(t)$ on the $SI$-plane whose tangent vector at some points
is different from the vector field at those points (in contrast with the
natural course). The distinction between $\bm{\gamma}$ and $\bm{\sigma}$ is
important, since the epidemic does not generally act on the population without
interventions and thus we observe an amended course. Two measures of
intervention effectiveness can be constructed based upon the difference of the
two courses. The first is based on the area between the positions under the
actual and the natural course during the time interval $(a,b)$. Thus, we
define

$L_{a,b}=\int\limits_{a}^{b}\frac{||\bm{\sigma}(t)-\bm{\gamma}(t)||}{||\bm{\sigma}(t)||}dt$
(16)

where $\bm{\sigma}$ is the natural and $\bm{\gamma}$ is the actual course.
Although the natural course is unknown in practice, setting it to the current
course and letting $\bm{\gamma}$ be the course under a hypothetical scenario
where an intervention measure is taken leads to a measure of effectiveness
quantification of this intervention. Applying the same idea for a different
kind of measure, allows one to compare the two interventions. Large values of
$L_{a,b}$ suggest large deviation of the course under the intervention from
the initial course. Thus, intervention A will be better than intervention B if
$L_{a,b}^{(A)}>L_{a,b}^{(B)}$. The second measure is obtained by comparing the
epidemic work between the natural and actual course during the time interval
$(a,b)$, where the epidemic work is defined as

$W_{a,b}=\int_{a}^{b}\left(\frac{dS}{dt}\right)^{2}+\left(\frac{dI}{dt}\right)^{2}dt$
(17)

i.e. the integrated squared speed

$v(t)=\sqrt{\left(\frac{dS}{dt}\right)^{2}+\left(\frac{dI}{dt}\right)^{2}}$
(18)

Lower speed is indicative of a better course, thus the larger the difference
between the actual course and the one suggested by the intervention, the
better. To this end, we define

$M_{a,b}=\frac{\displaystyle\int_{\bm{\sigma}}\mathbf{F}ds-\int_{\bm{\gamma}}\mathbf{F}dg}{\displaystyle\int_{\bm{\sigma}}\mathbf{F}ds}$
(19)

as the second intervention effectiveness measure. Again,
$M_{a,b}^{(A)}>M_{a,b}^{(B)}$ means intervention A is to be preferred over
intervention B.

Note that the simple SIR system (1) can lead to a conserved quantity
$Q(t)=S(t)+I(t)-\displaystyle\frac{1}{\lambda\tau}\log\big{(}S(t)\big{)}$,
which characterizes this particular system. Inspecting $Q(t)$ at the start of
the epidemic gives a measure of the deviation of the realised epidemic from
the SIR formulation. Under SIR, $Q_{0}(t)$ is constant for initial values
$(S_{0},I_{0},R_{0})$, so $(Q_{0}^{*}(t)-Q_{0}(t))^{2}$ (where $Q_{0}^{*}(t)$
is obtained by a different model) measures the deviation from the SIR for the
specific set of initial values $(S_{0},I_{0},R_{0})$ for a specific time $t$.
The absolute difference may also be used in this measure.

###  4.2 Estimation on the Discrete Time Domain

Since our models are defined in discrete time we develop here the discrete
analogues and illustrate these concepts on the daily Covid-19 data from
Greece. The SEIR with vaccination and demography model proposed in this
article provides estimates for the posterior medians of the susceptible and
infectious individuals, which we plot in Figure 5 on the $SI$-plane as
proportions of the total population. The red dots plot the estimated epidemic
course $\hat{\bm{\gamma}}_{t}=(\hat{S}_{t},\hat{I}_{t})$, is initiated from at
the bottom right corner
$\hat{\bm{\gamma}}_{0}=(N_{GR}-\hat{C}_{1},\hat{C}_{1})$. We include a top
axis whose ticks correspond to each $(\hat{S}_{t},\hat{I}_{t})$ point
facilitating for easy reading of the points of high and low speed (with
respect to changes in $S$) given in (18), which now reads

$v_{t}=\sqrt{(S_{t+1}-S_{t})^{2}+(I_{t+1}-I_{t})^{2}}$ (20)

while the epidemic work given in (17) can be obtained by

$W_{a,b}=\sum\limits_{t=a}^{b-1}(S_{t+1}-S_{t})^{2}+(I_{t+1}-I_{t})^{2}$ (21)

At the same Figure we see that the epidemic at the end of 2021 was at the
point $(S,I)=(0.4388,0.0039)$, so approximately 4746069 of the Greek
population escaped infection, while 421234 was infectious at that time.

Figure 5: The SEIR with vaccination and demography model estimates on the
$SI$-plane.

The intervention effectiveness for the time interval $(a,b)$ may be measured
using

$L_{a,b}=\sum\limits_{t=a}^{b}\sqrt{\frac{(S_{t}^{(n)}-S_{t}^{(a)})^{2}+(I_{t}^{(n)}-I_{t}^{(a)})^{2}}{S_{t}^{(n)2}+I_{t}^{(n)2}}}$
(22)

and

$M_{a,b}=\frac{W_{a,b}^{(n)}-W_{a,b}^{(a)}}{W_{a,b}^{(n)}}$ (23)

where the superscripts $(n)$ and $(a)$ indicate natural and actual course
respectively (or courses under the current and a new model).

Figure 6: Left: $Q_{t}$ for the $(S,I)$ values of Greece. There is a downward
pattern after approximately 224 days. Right: Zoomed in picture of the first
225 days for better visualization of the random fluctuations of $Q_{t}$,
albeit approximately constant in that interval. Red line is used for the
ergodic mean of the points.

We use the estimated $(S,I)$ quantities for Greece to examine the SIR
assumption by checking the quantity $Q_{t}=S_{t}+I_{t}-\log
S_{t}/(\lambda_{t}\tau)$. If the SIR assumption is correct, then $Q_{t}$ has
to be roughly constant. In the left plot of Figure 6 a downward trend is
apparent, invalidating the SIR assumption beyond the first 224 days. A zoomed
version for this time interval is depicted in the right plot of the same
Figure. Due to the random fluctuations in $S$ and $I$ we also plot the ergodic
mean (in red). The $Q_{t}$ values are around 1.000349 and this approach may
serve as a simple goodness of fit test for the SIR formulation.

##  5 Discussion

We found that the first two years of the Covid-19 epidemic may reasonably be
described by SEIR-type models which include information on demographic births
and deaths as well as vaccination. Our inference is based upon the estimated
total cases and the external validity of the results is inspected through the
comparison with the REACT-2 study in UK. The full SEIR model with vaccination
and demography effects is tested removing some of its components in turn and
comparisons are made using information criteria, the marginal likelihood and
Bayes factors. The proposed framework is developed based on publicly available
data which are the focus of this work. Estimation is based on the NUTS variant
of HMC since maximization of the un-normalized posterior via Simulated
Annealing appeared unstable. The commonly used Negative Binomial was assessed
and the Exponential-based alternative mixture offered a statistically robust
and computationally efficient alternative model. The predictive ability of the
mobility data is examined within a dimension reduction framework.

A potential weakness of the proposed model is the reliance on an IFR estimate
which affects some of the model parameters. This is an issue because we
focused upon publicly available data but our framework can accommodate
additional evidence, like seroprevalence surveys, in a straightforward manner.
Such additional data would give direct evidence on the IFR but are not
typically available for many countries and therefore we selected not to
utilise them. The computational burden for training our models poses limits to
certain short-term prediction exercises. However, one-week-ahead predictions
are certainly feasible and often the most realistic, without imposing strong
assumptions on the population behaviour. We selected not to use case data when
fitting the main model as we wanted to work with the estimated total number of
cases. One may add these data if interested, say, upon short term prediction,
but we selected to use case data post-hoc in order to estimate quantities like
the proportion observed.

The dynamics of the proposed epidemic models are examined on the phase plane
of the susceptible and infectious individuals, offering estimates of
effectiveness for control measures such as the non-pharmaceutical
interventions adopted in the early phase of the pandemic. A goodness-of-fit
type of assessment for the SIR assumption is proposed based on the time-
invariance of the $Q_{t}$ quantity. Further work based on variants of $Q_{t}$
is required and this is the subject of future research.

## References

  * (1) https://www.statista.com/statistics/270000/age-distribution-in-the-united-states/. Accessed in 2023.
  * (2) Apple (2020).  Apple mobility trends reports.
  * Cereda et al. (2020) Cereda, D., Tirani, M., Rovida, F., Demicheli, V., Ajelli, M., Poletti, P., …, Merler S. (2020)  The early phase of the Covid-19 outbreak in Lombardy, Italy.  arXiv preprint arXiv:2003.09320.
  * (4) https://www.cdc.gov/coronavirus/2019-ncov/hcp/planning-scenarios.html. Accessed in 2021.
  * (5) https://www.ethnicity-facts-figures.service.gov.uk/uk-population-by-ethnicity/demographics/age-groups/latest. Accessed in 2022.
  * Dong et al. (2020) Dong, E., Du, H. & Gardner, L. (2020).  An interactive web-based dashboard to track Covid-19 in real time.  The Lancet infectious diseases 20(5), 533–534.
  * Flaxman et al. (2020) Flaxman, S., Mishra, S., Gandy, A., Unwin, H. J. T., Mellan, T. A., Coupland, H, …, Bhatt S. (2020).  Estimating the effects of non-pharmaceutical interventions on Covid-19 in Europe.  Nature 584(7820), 257–261.
  * (8) Google (2020).  Google covid-19 community mobility reports.
  * He et al. (2020) He, X., Lau, E. H., Wu, P., Deng, X., Wang, J., Hao, X., …, Leung, G. M. (2020)  Temporal dynamics in viral shedding and transmissibility of Covid-19.  Nature Medicine 26(5), 672–675.
  * Hoffman et al. (2014) Hoffman, M. D. & Gelman, A. (2014).  The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research 15 15(1), 1593–1623.
  * (11) https://www.statistics.gr/el/statistics/-/publication/SAM03/-. Accessed in 2022.
  * Lauer et al. (2020) Lauer, S. A., Grantz, K. H., Bi, Q., Jones, F. K., Zheng, Q., Meredith, H. R., …, Lessler J. (2020).  The incubation period of Coronavirus disease 2019 (Covid-19) from publicly reported confirmed cases: estimation and application.  Annals of internal medicine 172(9), 577–582.
  * Mathieu et al. (2021) Mathieu, E., Ritchie, H., Ortiz-Ospina, E., Roser, M., Hasell, J., Appel, C, …, Rodés-Guirao, L. (2021).  A global database of Covid-19 vaccinations.  Nature human behaviour 5(7), 947–953.
  * Paul and Lorin (2021) Paul, S. & Lorin, E. (2021).  Estimation of Covid-19 recovery and decease periods in Canada using delay model.  Scientific Reports 11(1), 1–15.
  * (15) https://coronavirus.data.gov.uk/. Accessed in 2022.
  * (16) https://github.com/Sandbird/covid19-Greece. Accessed in 2022.
  * (17) Stan Development Team (2022).  RStan: the R interface to Stan.  R package version 2.21.5.
  * (18) https://www.statista.com/statistics/263762/total-population-of-the-united-states/. Accessed in 2023.
  * (19) https://www.worldometers.info/world-population/uk-population/. Accessed in 2022.
  * Ward et al. (2021) Ward, H., Atchison, C., Whitaker, M., Ainslie, K. E., Elliott, J., Okell, L., …, Elliott, P. (2021).  SARS-CoV-2 antibody prevalence in England following the first peak of the pandemic.  Nature communications 12(1), 1–8.


# Leveraging a Joint of Phenotypic and Genetic Features on Cancer Patient
Subgrouping

DODavid Oniani  CWChen Wang  YZYiqing Zhao  WAAndrew Wen  HLHongfang Liu
FSFeichen Shen  Kern Center for the Science of Health Care Delivery, Mayo
Clinic, Rochester, US  Division of Biomedical Statistics and Informatics, Mayo
Clinic, Rochester, US  Division of Digital Health Sciences, Mayo Clinic,
Rochester US

###### Abstract

Background Cancer is responsible for millions of deaths worldwide every year.
Although significant progress has been achieved in cancer medicine, many
issues remain to be addressed for improving cancer therapy. Appropriate cancer
patient stratification is the prerequisite for selecting appropriate treatment
plan, as cancer patients are of known heterogeneous genetic make-ups and
phenotypic differences. In this study, built upon deep phenotypic
characterizations extractable from Mayo Clinic electronic health records
(EHRs) and genetic test reports for a collection of cancer patients, we
developed a system leveraging a joint of phenotypic and genetic features for
cancer patient subgrouping.

Methods The workflow is roughly divided into three parts: feature
preprocessing, cancer patient classification, and cancer patient clustering
based. In feature preprocessing step, we performed filtering, retaining the
most relevant features. In cancer patient classification, we utilized joint
categorical features to build a patient-feature matrix and applied nine
different machine learning models, Random Forests (RF), Decision Tree (DT),
Support Vector Machine (SVM), Naive Bayes (NB), Logistic Regression (LR),
Multilayer Perceptron (MLP), Gradient Boosting (GB), Convolutional Neural
Network (CNN), and Feedforward Neural Network (FNN), for classification
purposes. Finally, in the cancer patient clustering step, we leveraged joint
embeddings features and patient-feature associations to build an undirected
feature graph and then trained the cancer feature node embeddings.

Results In cancer patient classification, the optimal performance of F1 0.9
was achieved by utilizing joint features. For cancer feature graph embeddings,
average operation for edge embeddings yielded the highest ROC score of 0.71,
followed by Hadamard. Regarding cancer patient clustering, CNN in combination
with fine-tuned cancer feature graph embeddings (CFEmb+) showed the best
purity score of 0.48, followed by UMLS/Mut2vec embeddings (UMEmb). In the
qualitative evaluation, we observed some generated clusters had known
relationships between diseases, phenotypes and genetic mutations.

Conclusions We utilized Mayo Clinic EHRs and genetic test reports in order to
design a system for cancer patient subgrouping. We discovered the potential
usage of applying machine learning algorithms over joint phenotypic and
genetic features for a better cancer type characterization.

Cancer,

Patient Subgrouping,

Joint Feature Analysis,

Machine Learning,

###### keywords:

Research

##  1 Introduction

Cancer is a type of disease in which abnormal cells divide without control and
can invade nearby tissues and other parts of the body via the blood and lymph
systems [1, 2]. Nowadays, cancer is responsible for millions of deaths
worldwide every year. Although significant progress has been achieved in
cancer medicine, many issues remain to be addressed for improving cancer
therapy. The goal of Precision Oncology is to enable oncologist practitioners
to make better clinical decisions by incorporating individual cancer patients’
genetic information and clinical characteristics [3]. Precision Oncology is
expected to improve selection of targeted therapies, tailored to individual
patients and ultimately improve cancer patients’ outcomes. However, patient
heterogeneity needs to be considered while providing cancer treatments. It is
important to take patient subgroups into consideration in which the mechanism
of diseases with the same group are more likely to be homogeneous [4].
Therefore, appropriate cancer patient subgrouping is the prerequisite for
selecting appropriate treatment plan, as cancer patients are of known
heterogeneous genetic make-ups and phenotypic differences.

Advances in cancer care have demonstrated the benefits of Individualized
Medicine research, and the premise of Precision Oncology is to refer similar
cancer outcomes to similar patients with similar features [5]. Phenotypic
characterizations are one of the key features to group patients by analyzing
their clinical observations [6]. Taking breast cancer as an example, specific
sub-phenotypes have been leveraged to make decision making for therapies
selection in practice [7]. Moreover, phenotypic heterogeneity amongs 7
subtypes of triple-negative breast cancer could be observed and analyzed [8].
With the advance of machine learning techniques, computational deep
phenotyping is widely used and described as the process to analyze phenotypic
abnormalities precisely and comprehensively, in order to receive a better
understanding of the natural history of a disease and its genotype-phenotype
associations. For example, the DeepPhe information model is able to extract
cancer-related phenotypes as categorical features by mining clinical
narratives or structured datasets [9]. Beeghly-Fadiel et al. have extended
DeepPhe to ovarian cancer [10]. Besides categorical phenotypic features,
embeddings representations of phenotypes are also significant to deep
phenotyping. HPO2Vec+ provides a way to learn node embeddings for the Human
Phenotype Ontology (HPO) and be able to represent patients with vectorized
representations for individual-level phenotypic similarity measurement [11,
12]. Maldonado et al. presented an approach that uses Generative Adversarial
Networks (GANs) [13] for learning UMLS embeddings and demonstrated the use
case on clinical prediction tasks [14].

In addition to phenotypic characterizations, features contained in genetic
test reports provide crucial evidence for the decision making in practice.
Genetic tests for targeted cancer therapy detect mutations in the DNA of
cancer cells. Being aware of whether the cancer has a particular mutation can
help guide the type of treatment. [15]. Specifically, pathogenic variants
demonstrate the genetic alteration that increases an individual’s
susceptibility or predisposition to a certain cancer type, and variant of
uncertain significance (VUS) indicates a genetic variant that has been
identified but without knowing too much details regarding significance to the
function of an organism. There exist some studies that use machine learning to
support variant interpretation in a clinical setting. For example, the LEAP
model classifies variants based on features collected from functional
predictions, splice predictions, clinical observation data and so on [16].
Some guidelines to define variant pathogenicity for somatic and VUS
classification were established [17, 18, 19]. Mut2Vec provides a way to
represent cancerous mutations with distributed embeddings, in order to support
cancer classification and drug sensitivity prediction [20].

Although phenotypic and genetic features have their own advantages in
facilitating computational individualized medicine and precision oncology,
incorporating longitudinal clinical phenotypes into mutation features for
cancer patient stratification is largely unexplored. In this study, built upon
a cancer patient cohort with both deep phenothypic characterizations
extractable from Mayo Clinic electronic health records (EHRs) and well-
documented genetic test reports, we designed an approach leveraging a joint of
phenotypic and genetic features for cancer patient subgrouping.

In the following, we first introduce materials and methods used in this study.
We then describe experimental strategies followed by Results. We conclude our
study with discussion and future work.

##  2 Materials

We collected a cancer patient cohort (n=794) with accessible Mayo Clinic EHRs
and foundation genetic test reports. Seven types of cancer are included in
this cohort: lung cancer (n=223), prostate cancer (n=66), breast cancer
(n=53), ovarian cancer (n=91), pancreas cancer (n=104), colon/rectum cancer
(n=140), and liver cancer (n=107).

###  2.1 Extract phenotypic features from the Mayo Clinic EHRs

Given patient IDs retrieved from the cohort, we first applied an existing Mayo
Clinic NLP-as-a-service implementation [21] to extract all the diagnostic-
related terms from clinical notes and radiology report, normalized by the
Unified Medical Language System (UMLS) [22]. We then used a previously
developed Human Phenotype Ontology (HPO)-based annotation pipeline [23] to
further map the UMLS terms to more specific HPO normalized phenotypic terms
[24]. Specifically, for each cancer patient, we collected their phentoypic
information within the last twelve months of the first mention of cancer
diagnoses.

###  2.2 Collect genetic feature from the FoundationOne CDx genetic test
reports

FoundationOne CDx is a next-generation sequencing (NGS) based assay that
identifies genomic findings within hundreds of cancer-related genes. The
foundationOne genetic test report is designed to provide physicians with
clinically actionable information. The report for each patient is represented
in the PDF format. We first used a python library named pdfminer to parse all
the PDF report into text format, for the convenience of data processing. We
then parsed out genes alteration and VUS sections for each patient. In
particular, we treated mutations mined from the genes alteration section as
pathogenic mutation, and we named mutations mined from the VUS section as VUS
mutation. We used multi-hot encoding to represent genetic features for each
patient, where 1 indicates the existence of pathogenic mutation or VUS, and 0
otherwise.

##  3 Methods

The workflow of this study is roughly divided into three parts. As shown in
Figure Figures, first of all, we preprocessed the extracted phenotypic and
genetic features to remove noises. In addition, we adopted supervised learning
strategy over the patients for cancer patient classification. Furthermore, we
constructed a cancer feature graph based on associations between patients and
different features, and we trained node embeddings based on the feature graph
to get the vectorized feature representation. We then applied unsupervised
learning approaches over the patient-level embeddings for patient clustering.

###  3.1 Feature preprocessing

For phenotypes, we first removed some general phenotypic terms mined from
EHRs, including cyst, pain, carcinoma, neoplasm, symptoms, disease, minor
(disease), and sarcoma - cateogry (morphologic abnormality). We then filtered
out phenotypic terms with low frequency across all the patients, in order to
largely avoid noises. By heuristics, we set the frequency threshold for
selecting phenotypes as 20. Similarly, for genetic features, we only kept the
ones that appear more than 10 times across all the patients.

###  3.2 Cancer patient classification based on joint categorical features

We built a $M*N$ patient-feature matrix with multi-hot encoding, where $M$
indicates the number of patients and $N$ denotes the total number of features,
including 239 phenotypic features, 328 genetic features, and 1 demographic
feature for gender. Given patients’ cancer diagnosis as labels, we treated
cancer patient classification as a multi-class classification task and applied
different supervised learning classifiers over the patient-feature matrix to
categorize patients based on the diagnosed cancer types. Here, we used Random
Forests (RF), Decision Tree (DT), Support Vector Machine (SVM), Naive Bayes
(NB), Logistic Regression (LR), Multilayer Perceptron (MLP), Gradient Boosting
(GB), Convolutional Neural Network (CNN), and Feedforward Neural Network
(FNN).

###  3.3 Cancer patient clustering based on joint embeddings features

####  3.3.1 Construction of the cancer feature graph

As shown in Figure Figures, associations between the patients and the features
were captured to build the feature graph. Specifically, we made a link between
phenotypic feature and genetic feature if both of them are associated with the
same patient. As a result, the feature graph can be defined as $G=(N,L)$,
where $N$ represents the set of phenotypic/genetic nodes and $L$ represents
the set of links/edges between these nodes. We viewed graph $G$ as an
undirected graph that bridges mutations and phenotypes based on information
collected from heterogeneous cancer patients in clinical practice.

####  3.3.2 Train cancer feature node embeddings

We proceed by training node embeddings for feature node contained in the
cancer feature graph $G$.

Due to its demonstrated high performance on link prediction, we applied
node2vec model [25]. It implements a second-order random walk over the
topological structure of the graph. In every walk, there are three types of
nodes involved. Specifically, source, intermediate, and target nodes which in
these case all are phenotypic/genetic nodes.

In graph $G$, given any source entity node as $E_{s}$, target entity node as
$E_{t}$, we denote the intermediate entity that exists on the path between
$E_{s}$ and $E_{t}$ as $E_{st}$. Then the distribution of entity $E_{t}$ with
a fixed length of random walk can be represented as follows:

$P(E_{t}\mid E_{i})=\systeme*{\dfrac{\pi(E_{i},E_{t})}{Z}\ if\ (E_{i}\,E_{t})\
is\ a\ link,0\ otherwise}$ (1)

where $Z$ is the normalization constant and $\pi(E_{i},E_{t})$ is the
transition probability between the phenotypic/genetic nodes $E_{i}$ and
$E_{t}$.

We denote the weight over the link $(E_{i},E_{t})$ as $w_{(}E_{i},E_{t})=1$.
Then $\pi(E_{i},E_{t})$ can be calculated as follows:

$\pi(E_{i},E_{t})=\alpha(E_{s},E_{t})\cdot w(E_{i},E_{t})$ (2)

where $\alpha(E_{s},E_{t})$ is a search bias term and as defined in node2vec,
can be computed:

$\alpha(N_{s},N_{t})=\systeme*{\dfrac{1}{q}\ if\ D=2,1\ if\ D=1,\dfrac{1}{p}\
if\ D=0}$ (3)

where $D$ is the shortest distance between the nodes $E_{s}$ and $E_{t}$ and
$p$ and $q$ are hyperparameters that balance between the breadth-first search
and depth-first search searching strategies for both local and global
optimization.

After learning the sampled network data using random walk, we then leveraged
the Skip-Gram model to train entity representations on the sampled data. For
each entity node$E_{s}\in N$and all its sampled neighbors $SN(E_{s})$, the
loss function for entity representation learning could be described as shown
in Equation (4). This loss function via Stochastic gradient descent. The
prediction distribution was normalized using softmax function.

$\underset{f}{\max}=\sum_{E_{s}\in N}\log P(SN(E_{s})\mid f(E_{s}))$ (4)

As a result, we have mapped the bipartite cancer feature graph from graph
space to embeddings space. Each node is represented by a trained vector and
strength of association for each pair of node is measurable.

####  3.3.3 Patient-level embeddings representation

For any patient $p$, the final patient level embeddings with respect to
phenotype or genetic features will be represented by using average embeddings
(AE):

$AE(p)=\dfrac{\sum_{me_{p}\in pv_{p}}^{|pv_{p}|}f_{c}(me_{p})}{|pv_{p}|}$ (5)

where $f_{c}$ a feature matrix to map any phenotype or mutation $me_{p}$
contained in patient vector $pv_{p}$ to its corresponding embeddings, and $c$
indicates which graph embeddings to be used (either our trained embeddings or
existing ones). Given the initial patient embeddings with respect to phenotype
as $R_{p}$, patient embeddings with respect to genetic features as $R_{g}$,
and demographic multi-hot vector as $R_{d}$, the final patient vector $R_{f}$
will be represented by concatenating three vectors together as
$[R_{p},R_{g},R_{d}]$

####  3.3.4 Patient clustering

We further fine-tuned the initial patient-level embeddings using the cancer
classification task via a categorical cross-entropy loss. We then adopted the
t-distributed stochastic neighbor embedding (t-SNE) algorithm to render the
refined patient embeddings for all patients into a 2D space [26]. We first
visualized patient subgrouping intuitively with t-SNE. As t-SNE does not
perform clustering in and of itself, but instead renders each patient
embedding into a $(x,y)$ coordinate. As such, additional clustering approach
is required to re-categorize these points into discrete clusters. K-means
clustering [27] was therefore applied over patients to further partition
different patient groups into distinct clusters.

##  4 Evaluation

###  4.1 Evaluate supervised learning performance

We used 70% data for training purpose and the rest for testing. We conducted a
10-fold cross-validation approach over the training data for choosing the
optimal hyper-parameters for each classifier. To handle data imbalanced
problem, we oversampled the training data. In order to demonstrate the
advantage of combining both phenotypic and genetic features for cancer patient
classification, we compared the cancer patient classification performance by
using phenotypic feature only, genetic feature only, and the joint of two type
of features. We used precision, recall, and F1-score to evaluate the overall
classification performance and which for each specific cancer.

###  4.2 Evaluate the construction of node embeddings

We leveraged a link prediction task to train node embeddings for the cancer
feature graph. We used 60% of edges as training data, 10% edges as validation
set, and the rest is used for testing purpose. Similar to what reported in the
node2vec study, we tested edge embeddings with four operations, including
Hadamard, Average, L1, and L2 [25]. We combined RF with different edge
embeddings to train the optimal node embedding through the link prediction
task, and we reported ROC score for each experiment.

###  4.3 Evaluate unsupervised learning performance

We compared the trained cancer feature graph embeddings with other existing
embeddings on cancer patient clustering. Specifically, UMLS embeddings [14]
and Mut2vec [20] are two pre-trained embeddings for representing phenotype and
genetic mutation respectively. In this evaluation, according to Equation (5),
we used features contained in UMLS embeddings and Mut2vec to represent
patient-level embeddings, and then compared patient subgrouping performance
generated by six experimental settings: original cancer feature graph
embeddings (CFEmb), fine-tuned cancer feature graph embeddings with CNN
(CFEmb+ with CNN), fine-tuned cancer feature graph embeddings with FNN (CFEmb+
with FNN), original UMLS/Mut2vec embeddings (UMEmb), fine-tuned UMLS/Mut2vec
embeddings with CNN (UMEmb+ with CNN), and fine-tuned UMLS/Mut2vec embeddings
with FNN (UMEmb+ with FNN). Size of dimension of the pre-trained UMLS
embeddings and Mut2vec are 51 and 300 respectively. Since UMLS embeddings and
Mut2vec don’t encode gender as a feature, we used a 1-bit binary indicator to
represent gender here. Therefore, the total dimension size for patient-level
embeddings represented by UMLS embeddings and Mut2vec is 352. We visualized
the t-SNE plots and selected the optimal number of clusters (ranging from 2 to
20) by measuring the purity scores for K-means clustering across all the
experiments. We also conducted a qualitative analysis for clusters generated
by the experiment with the optimal purity score. We picked the embeddings with
the optimal purity score and analyzed the most frequent cancers, phenotypes,
and genetics involved in some selected clusters.

##  5 Results

###  5.1 Cancer patient classification

As shown in Table 1, for CNN, FNN, LR, MLP, RF, GB, and SVM, performance using
joint features outpeformed the usage of phenotypic and genetic features alone
respectively. We observed that phenotypic features collected from EHRs are
able to better differentiate patients than genetic features collected from
genetic test report in this study. In particular, for DT and NB, using
phenotypic features alone outperformed using the combination of features. For
joint features, LR and MLP outperformed the rest of the models, with LR
showing marginal F1-score advantage over MLP. FNN, RF, GB, and SVM also
performed similarly with the performance slightly lower than both LR and MLP.
NB classifier showed the worst performance with the precision, recall, and F1
values of 0.48, 0.40, and 0.40 respectively. CNN showed the second worst
performance with joint multi-hot encoding features. Its precision value was
0.77, recall was 0.73, and F1 score was 0.74.

Table 1: Comparison of overall cancer patient classification performance between different features. The optimal precision, recall and F1 scores are in bold. Model | Joint features | Phenotypic features | Genetic features
---|---|---|---
| Precision | Recall | F1 Score | Precision | Recall | F1 Score | Precision | Recall | F1 Score
CNN | 0.77 | 0.73 | 0.74 | 0.75 | 0.74 | 0.73 | 0.24 | 0.37 | 0.25
FNN | 0.88 | 0.86 | 0.87 | 0.85 | 0.86 | 0.85 | 0.15 | 0.19 | 0.17
DT | 0.76 | 0.78 | 0.76 | 0.83 | 0.77 | 0.77 | 0.15 | 0.19 | 0.17
LR | 0.90 | 0.90 | 0.90 | 0.84 | 0.84 | 0.83 | 0.21 | 0.37 | 0.25
MLP | 0.89 | 0.88 | 0.88 | 0.83 | 0.82 | 0.82 | 0.19 | 0.36 | 0.23
NB | 0.48 | 0.40 | 0.40 | 0.52 | 0.54 | 0.51 | 0.24 | 0.29 | 0.25
RF | 0.87 | 0.87 | 0.86 | 0.85 | 0.85 | 0.83 | 0.19 | 0.25 | 0.21
GB | 0.89 | 0.85 | 0.86 | 0.86 | 0.85 | 0.85 | 0.19 | 0.24 | 0.21
SVM | 0.90 | 0.86 | 0.87 | 0.87 | 0.86 | 0.85 | 0.21 | 0.25 | 0.22

Table 2 is in accordance with the overall results shown in Table 1. The
average F1 score for the classification of ovarian cancer was the highest,
equal to 0.85. Liver cancer came next, with the average F1 score of 0.81.
Average F1 score for the breast cancer classification was 0.78. Prostate
cancer and pancreas cancer had the average F1 scores of 0.77 and 0.74
respectively. Lung cancer had the average F1 score of 0.73 and NB was not able
to make any correct classification for lung cancer patients. Colon/rectum
cancer had the worst F1 score across all models used with joint feature
representations, with the score of 0.72.

Table 2: Classification F1 score on each cancer type using joint features. The optimal F1 scores are in bold. Model | Lung | Prostate | Breast | Ovarian | Pancreas | Col./Rec. | Liver
---|---|---|---|---|---|---|---
CNN | 0.75 | 0.85 | 0.76 | 0.85 | 0.76 | 0.75 | 0.87
DNN | 0.67 | 0.68 | 0.45 | 0.81 | 0.60 | 0.55 | 0.68
DT | 0.67 | 0.62 | 0.83 | 0.81 | 0.75 | 0.70 | 0.83
LR | 0.89 | 0.92 | 0.88 | 0.93 | 0.88 | 0.88 | 0.91
MLP | 0.89 | 0.93 | 0.82 | 0.93 | 0.88 | 0.81 | 0.91
NB | 0.00 | 0.43 | 0.60 | 0.61 | 0.55 | 0.25 | 0.37
RF | 0.84 | 0.87 | 0.88 | 0.92 | 0.75 | 0.86 | 0.90
GB | 0.89 | 0.86 | 0.87 | 0.87 | 0.74 | 0.86 | 0.93
SVM | 0.94 | 0.81 | 0.89 | 0.93 | 0.77 | 0.84 | 0.91

###  5.2 Cancer feature graph embeddings

As shown in Figure Figures, link prediction performance made by RF with four
edge embeddings operations were compared. We found that Average edge
embeddings yielded the highest ROC score of 0.71, followed by Hadamard (area =
0.70). L1 and L2 had the same prediction performance of 0.63, respectively.
Therefore, we used RF with Average edge embeddings to train the embeddings and
the final hyper-parameters were selected as follows: return parameter $p$ is
0.5, in-out parameter $q$ is 0.5, window size is 10, number of walk is 5, walk
length is 10, dimension = 100.

###  5.3 Cancer patient clustering

We first conducted a qualitative analysis of the patient-level clustering. Six
t-SNE clustering visualizations were made as shown in Figure Figures.
Regarding the usage of our self-trained cancer feature graph embeddings, we
observed that the original CFEmb was not able to subgroup cancer patients
clearly based on the joint embeddings features. However, we found that fine-
tuned CFEMb with both CNN and FNN were able to better organize patients with
same cancer types. Particularly, there is no significant differences between
CFEmb+ with CNN and CFEmb+ with FNN could be observed from visualizations.
Intuitively, we found that refined features from both CFEmb+ can detect clear
clusters for patients with prostate cancer (blue), breast cancer (green),
ovarian cancer (black), and liver cancer (red). For the combination usage of
UMLS embeddings and Mut2vec, we observed that there exists no significant
differences on clustering cancer patients among the original UMEmb, UMEmb+
with CNN, and UMEmb+ with FNN. Compared to CFEmb+, UMEmb+ are less capable of
making cohesion clusters to represent each type of cancer. Moreover, since
CFEmb+ encoded gender as embeddings, we found that both CFEmb+ with CNN and
FNN were able to make a clear separation for patients with different genders.
For example, for CFEmb+ with CNN, all the male patients were grouped in the
left part of the visualization (e.g., patients with prostate cancer (blue))
while all the female patients were grouped in the right part, including
patients with breast cancer (green) and ovarian cancer (black). The same
separation could be observed in CFEmb+ with FNN. While neither UMEmb+ with CNN
nor UMEmb+ with FNN were capable of separating patients by gender.

In addition, we made a quantitative analysis for six embeddings using K-means.
Table 3 depicts the number of clusters along with the associated best purity
score selected for each embeddings. We observed that CFEmb+ with CNN
conbritued to the optimal purity score of 0.48 while $k=18$.

Table 3: Clustering performance over different patient-level embeddings. The optimal purity score is in bold. Methods | Number of clusters (k) | Purity score
---|---|---
CFEmb | 12 | 0.30
CFEmb+ with FNN | 19 | 0.39
CFEmb+ with CNN | 18 | 0.48
UMEmb | 18 | 0.41
UMEmb+ with FNN | 18 | 0.39
UMEmb+ with CNN | 17 | 0.39

We picked CFEmb+ with CNN as the optimal embeddings and further analyzed some
selected clusters. We observed that several generated clusters had known
relationships between diseases, phenotypes and genetic mutations. For example,
as shown in Table 4, cluster-8 was mostly enriched with prostate cancer (n=39)
and related disease phenotypes (prostatic disease, carcinoma in situ of
prostate), as well as TP53 pathogenic mutations known implicated in prostate
cancer progression [28]; in addition, pathogenic alterations of androgen
receptor (AR) is found associated with cluster-8, and has been clearly
established as driver genetic event leading to prostate cancer and considered
as target for various therapeutic options for treatment prostate cancer [29].
As two other examples, cluster-7 and cluster-3 were both enriched with
colon/rectum cancers (n=24 and 21 respectively) and closely related phenotypes
(e.g. colonic diseases); in addition, commonly mutated TP53 and KRAS genes
were found in both clusters, with KRAS alternation as an established
treatment-response marker and potential drug targets [30, 31]. Interestingly,
cluster-7 and cluster-3 lead to uniquely different VUS findings in BRCA2
(cluster-7) and TSC2 (cluster-3), raising the scientific possibilities that
these colorectal disease were from different disease etiologies and mutational
background. Although these observations require larger sample-size and
additional cohort for independent validations, data-driven discovery procedure
as such clearly pave a way for formulating novel hypothesis through joint
mining of cancer, phenotypes and genetic changes.

Table 4: Selected clusters with top most frequent cancer type, phenotypic features, and genetic features. Number in bracket indicates the frequency of entity. ”-PATH” indicates a pathogenic mutation, while ”-VUS” indicates a variant of uncertain significance. Cluster ID | Cancer type | Phenotypic features | Genetic features
---|---|---|---
Cluster-8 | Prostate cancer (39) | Prostatic diseases (44) | TP53-PATH (18)
Colon/Rectum cancer (3) | Carcinoma in situ of prostate (44) | PTEN-PATH (17)
Lung cancer (2) | Neoplasm of uncertain or unknown behavior of prostate (44) | TAF1-VUS (18)
Pancreas cancer (2) | Benign neoplasm of prostate (44) | AR-PATH (18)
Cluster-7 | Colon/Rectum cancer (24) | Neoplasm of uncertain or unknown behavior of colon (35) | TP53-PATH (20)
Lung cancer (7) | Colonic diseases (35) | KRAS-PATH (20)
Pancreas cancer (3) | Stage 0 carcinoma of colon (35) | ASXL1-VUS (18)
Liver cancer (1) | Benign neoplasm of liver (28) | BRCA2-VUS (7)
Cluster-3 | Colon/Rectum cancer (21) | Neoplasm of uncertain or unknown behavior of colon (34) | TP53-PATH (33)
Lung cancer (6) | Colonic diseases (34) | ASXL1-VUS (18)
Ovarian cancer (4) | Stage 0 carcinoma of colon (34) | KRAS-PATH (15)
Pancreas cancer (4) | Urologic diseases (30) | TSC2-VUS (13)

##  6 Discussion and Future Work

For cancer patient classification with joint categorical features, we found
that CNN was not able to outperform FNN and other conventional machine
learning algorithms. The reason behind this is that the size of cancer patient
cohort is limited and the CNN models trained on the small dataset usually
yield poor generalization and horrible performance [32]. In addition, we
observed that NB with phenotypic features outperformed which with joint
features, despite the low F1 score of 0.51. Specifically, we used the
Gaussian-based NB in evaluation and we measured average feature importance
score using permutation feature importance measurement. We found that NB
achieved an average score of 0.0008 and 0.0023 by using joint features and
phenotypic features respectivley, indicating that homogeneous features were
more suitable for NB in this evaluation.

Regarding cancer feature graph generation, there is a debating between whether
to include patient as a graph node or not. We decided to not include patients
as nodes but only use them for bridging heterogeneous phenotypic and genetic
features, and here are rationales to implement it this way: 1) our plan is to
build a generic embeddings regarding important associations between phenotypes
and mutations. The constructed embeddings incorporates rich semantic
information contained in both Mayo Clinic EHRs and FoundationOne CDx genetic
test reports, which could be used as a pre-trained cancer feature embeddings
for othe datasets and studies independently; 2) Since we adopted the 2nd order
random walk within node2vec to train embeddings, adding patients as nodes in
between phenotypic and genetic features will reduce the chance of reachability
between phenotype and mutation nodes, thus weaken the power of association
prediction/detection amongs features.

For cancer patient clustering, in Table 4, we found that some most frequent
phenotypic features have coarse granularity. For example, in Cluster-8, all
the top 4 phenotypes contain keyword prostate or prostatic, that explains why
prostate cancer is the most frequent cancer type in this cluster. In the
future, we will explore more over phenotypic features with finer granularity,
in order to increase the model prediction power while dealing with phenotypic
features without any keyword mentions. To achieve this, we will focus on the
semantic relationship amongst synonym, hyponymy, and hypernymy within the HPO.
Our prior designed HPO2Vec+ framework [11] will also be leveraged to capture
strong semantic relationships for this purpose.

##  7 Conclusion

We utilized Mayo Clinic electronic health records (EHRs) and genetic test
reports in order to design a system for cancer patient subgrouping. For this
purpose, we developed a three-step process comprised of feature preprocessing,
cancer patient classification, and cancer patient clustering. We discovered
that the usage of the joint phenotypic and genetic features showed a
significant performance increase over the either set of features used alone.
The usage of genetic features alone showed the worst performance. In this
study, we discovered the potential usage of applying machine learning
algorithms over joint phenotypic and genetic features for a better cancer type
characterization.

## Acknowledgements

This work has been supported by the Gerstner Family Foundation, Center for
Individualized Medicine (CIM) of Mayo Clinic, and Genentech Research Fund in
Individualized Medicine.

## Abbreviations

EHRs: Electronic Health Records; RF: Random Forests; DT: Decision Tree; SVM:
Support Vector Machine; NB: Naive Bayes; LR: Logistic Regression; MLP:
Multilayer Perceptron; GB: Gradient Boosting; CNN: Convolutional Neural
Network; FNN: FeedforwardNeural Network; CFEmb: Cancer Feature Graph
Embeddings; CFEmb+: Fine-tuned Cancer Feature Graph Embeddings (CFEMb+);
UMEmb: UMLS/Mut2vec embeddings; UMEmb+: Fine-tuned UMLS/Mut2vec embeddings.

## Availability of data and materials

The cancer network file is stored in .txt format and the cancer network
embeddings file is saved in .emb format. Location of the data is at:
https://github.com/shenfc/Joint-Learning. Data files could be downloaded
freely from GitHub.

## Ethics approval and consent to participate

This study used existing data to conduct a retrospective study. The study and
a waiver of informed consent were approved by Mayo Clinic Institutional Review
Board in accordance with 45 CFR 46.116 (Approval #17–003030, #20-001474)All
methods were carried out in accordance with relevant guidelines and
regulations.

## Competing interests

The authors declare that they have no competing interests.

## Consent for publication

The author(s) declare(s) that the manuscript does not contain any individual
person’s data. So this paper requires no consent to publish.

## Authors’ contributions

All co-authors are justifiably credited with authorship, according to the
authorship criteria. Final approval is given by each co-author. In detail: DO
contributed to the system implementation, data analysis, result
interpretation, and manuscript preparation. CW contributed to system design,
result interpretation, and manuscript preparation. YZ contributed to data
extraction and data cleaning. AW is responsible for building data engineering
pipeline and extracting phenotypic data from electronic health records. HL
contributed to system design and evaluation. FS conceptualized and led the
study, contributed to system design and implementation, data analysis, results
interpretation, and manuscript preparation.

## References

  * [1] Yuan, Y., Shi, Y., Li, C., Kim, J., Cai, W., Han, Z., Feng, D.D.: Deepgene: an advanced cancer type classifier based on deep learning and somatic point mutations. BMC bioinformatics 17(17), 476 (2016)
  * [2] Feuerstein, M.: Defining cancer survivorship. Journal of cancer survivorship: research and practice 1(1), 5–7 (2007)
  * [3] Schwartzberg, L., Kim, E.S., Liu, D., Schrag, D.: Precision oncology: who, how, what, when, and when not? American Society of Clinical Oncology Educational Book 37, 160–169 (2017)
  * [4] on Neuroscience, F., Disorders, N.S., et al.: Improving and accelerating therapeutic development for nervous system disorders: Workshop summary (2014)
  * [5] Parimbelli, E., Marini, S., Sacchi, L., Bellazzi, R.: Patient similarity for precision medicine: A systematic review. Journal of biomedical informatics 83, 87–96 (2018)
  * [6] Seligson, N.D., Warner, J.L., Dalton, W.S., Martin, D., Miller, R.S., Patt, D., Kehl, K.L., Palchuk, M.B., Alterovitz, G., Wiley, L.K., et al.: Recommendations for patient similarity classes: results of the amia 2019 workshop on defining patient similarity. Journal of the American Medical Informatics Association (2020)
  * [7] Lehmann, B.D., Bauer, J.A., Chen, X., Sanders, M.E., Chakravarthy, A.B., Shyr, Y., Pietenpol, J.A., et al.: Identification of human triple-negative breast cancer subtypes and preclinical models for selection of targeted therapies. The Journal of clinical investigation 121(7), 2750–2767 (2011)
  * [8] Masuda, H., Baggerly, K.A., Wang, Y., Zhang, Y., Gonzalez-Angulo, A.M., Meric-Bernstam, F., Valero, V., Lehmann, B.D., Pietenpol, J.A., Hortobagyi, G.N., et al.: Differential response to neoadjuvant chemotherapy among 7 triple-negative breast cancer molecular subtypes. Clinical cancer research 19(19), 5533–5540 (2013)
  * [9] Savova, G.K., Tseytlin, E., Finan, S., Castine, M., Miller, T., Medvedeva, O., Harris, D., Hochheiser, H., Lin, C., Chavan, G., et al.: Deepphe: a natural language processing system for extracting cancer phenotypes from clinical records. Cancer research 77(21), 115–118 (2017)
  * [10] Beeghly-Fadiel, A., Warner, J.L., Finan, S., Masanz, J., Hochheiser, H., Savova, G.: Deep phenotype extraction to facilitate cancer research: Extending DeepPhe to ovarian cancer. AACR (2019)
  * [11] Shen, F., Peng, S., Fan, Y., Wen, A., Liu, S., Wang, Y., Wang, L., Liu, H.: Hpo2vec+: Leveraging heterogeneous knowledge resources to enrich node embeddings for the human phenotype ontology. Journal of biomedical informatics 96, 103246 (2019)
  * [12] Shen, F., Liu, S., Wang, Y., Wang, L., Wen, A., Limper, A.H., Liu, H.: Constructing node embeddings for human phenotype ontology to assist phenotypic similarity measurement. In: 2018 IEEE International Conference on Healthcare Informatics Workshop (ICHI-W), pp. 29–33 (2018). IEEE
  * [13] Goodfellow, I.J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial networks. arXiv preprint arXiv:1406.2661 (2017)
  * [14] Maldonado, R., Yetisgen, M., Harabagiu, S.M.: Adversarial learning of knowledge embeddings for the unified medical language system. AMIA Summits on Translational Science Proceedings 2019, 543 (2019)
  * [15] Sawyers, C.: Targeted cancer therapy. Nature 432(7015), 294–297 (2004)
  * [16] Lai, C., Zimmer, A.D., O’Connor, R., Kim, S., Chan, R., van den Akker, J., Zhou, A.Y., Topper, S., Mishne, G.: Leap: Using machine learning to support variant classification in a clinical setting. Human Mutation 41(6), 1079–1090 (2020)
  * [17] Federici, G., Soddu, S.: Variants of uncertain significance in the era of high-throughput genome sequencing: a lesson from breast and ovary cancers. Journal of Experimental & Clinical Cancer Research 39(1), 1–12 (2020)
  * [18] Richards, S., Aziz, N., Bale, S., Bick, D., Das, S., Gastier-Foster, J., Grody, W.W., Hegde, M., Lyon, E., Spector, E., et al.: Standards and guidelines for the interpretation of sequence variants: a joint consensus recommendation of the american college of medical genetics and genomics and the association for molecular pathology. Genetics in medicine 17(5), 405–423 (2015)
  * [19] Rehm, H.L., Bale, S.J., Bayrak-Toydemir, P., Berg, J.S., Brown, K.K., Deignan, J.L., Friez, M.J., Funke, B.H., Hegde, M.R., Lyon, E.: Acmg clinical laboratory standards for next-generation sequencing. Genetics in medicine 15(9), 733–747 (2013)
  * [20] Kim, S., Lee, H., Kim, K., Kang, J.: Mut2vec: distributed representation of cancerous mutations. BMC medical genomics 11(2), 33 (2018)
  * [21] Wen, A., Fu, S., Moon, S., El Wazir, M., Rosenbaum, A., Kaggal, V.C., Liu, S., Sohn, S., Liu, H., Fan, J.: Desiderata for delivering nlp to accelerate healthcare ai advancement and a mayo clinic nlp-as-a-service implementation. NPJ digital medicine 2(1), 1–7 (2019)
  * [22] Bodenreider, O.: The unified medical language system (umls): integrating biomedical terminology. Nucleic acids research 32(suppl_1), 267–270 (2004)
  * [23] Shen, F., Wang, L., Liu, H.: Phenotypic analysis of clinical narratives using human phenotype ontology. Studies in health technology and informatics 245, 581 (2017)
  * [24] Robinson, P.N., Köhler, S., Bauer, S., Seelow, D., Horn, D., Mundlos, S.: The human phenotype ontology: a tool for annotating and analyzing human hereditary disease. The American Journal of Human Genetics 83(5), 610–615 (2008)
  * [25] Grover, A., Leskovec, J.: node2vec: Scalable feature learning for networks. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 855–864 (2016)
  * [26] Maaten, L.v.d., Hinton, G.: Visualizing data using t-sne. Journal of machine learning research 9(Nov), 2579–2605 (2008)
  * [27] MacQueen, J., et al.: Some methods for classification and analysis of multivariate observations. In: Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, vol. 1, pp. 281–297 (1967). Oakland, CA, USA
  * [28] Ecke, T.H., Schlechte, H.H., Schiemenz, K., Sachs, M.D., Lenk, S.V., Rudolph, B.D., Loening, S.A.: Tp53 gene mutations in prostate cancer progression. Anticancer research 30(5), 1579–1586 (2010)
  * [29] Fujita, K., Nonomura, N.: Role of androgen receptor in prostate cancer: a review. The world journal of men’s health 37(3), 288–295 (2019)
  * [30] Siddiqui, A.D., Piperdi, B.: Kras mutation in colon cancer: a marker of resistance to egfr-i therapy. Annals of surgical oncology 17(4), 1168–1176 (2010)
  * [31] McCormick, F.: KRAS as a therapeutic target. AACR (2015)
  * [32] Zhang, G., Kato, J., Wang, Y., Mase, K.: How to initialize the cnn for small datasets: Extracting discriminative filters from pre-trained model. In: 2015 3rd IAPR Asian Conference on Pattern Recognition (ACPR), pp. 479–483 (2015). IEEE

## Figures

Figure 1. The overall workflow captionThe overall workflow

Figure 2. Building cancer feature graph captionBuilding cancer feature graph

Figure 3. ROC for link prediction captionROC for link prediction

Figure 4. t-SNE clustering for different embeddings captiont-SNE clustering
for different embeddings