# A new harmonium for pattern recognition in survival data

Hylke C. Donker<EMAIL_ADDRESS>Harry J. M. Groen  Department of Pulmonary
Diseases, University Medical Center Groningen, Hanzeplein 1, p.o. box 30.001,
9700 RB Groningen, the Netherlands

###### Abstract

Survival analysis concerns the study of timeline data where the event of
interest may remain unobserved (i.e., censored). Studies commonly record more
than one type of event, but conventional survival techniques focus on a single
event type. We set out to integrate both multiple independently censored time-
to-event variables as well as missing observations. An energy-based approach
is taken with a bi-partite structure between latent and visible states,
commonly known as harmoniums (or restricted Boltzmann machines). The present
harmonium is shown, both theoretically and experimentally, to capture non-
linear patterns between distinct time recordings. We illustrate on real world
data that, for a single time-to-event variable, our model is on par with
established methods. In addition, we demonstrate that discriminative
predictions improve by leveraging an extra time-to-event variable. In
conclusion, multiple time-to-event variables can be successfully captured
within the harmonium paradigm.

###### keywords:

survival analysis , machine learning , harmonium , restricted Boltzmann
machine

###### MSC:

[2020] 62N01 , 62H30 , 92B20

##  1 Introduction

Survival analysis considers the timing of dichotomous events, and can be used
to analyse, for example, time to death, breakdown of a machine, graduating
university or worsening of the disease. What distinguishes time-to-event
measurements from other random variables is that they are partially observed.
That is to say, not all events have occurred during the interval in which they
were monitored (we have lower bounds) because, e.g., subjects may be lost
during follow-up, or the experiment may be too short to observe all the
events. These incomplete measurements, where the event time remains unknown,
are said to be censored. Nevertheless, the time interval in which subjects
were observed prior to dropout still provides valuable information. In fact,
failure to take into account censoring leads to serious underestimation of the
survival, as has been repeatedly emphasised [1, 2, 3, 4].

A wide range of statistical tools have been developed to deal with censored
data, which rely on modelling the survival distribution
$S(t)=\int_{t}^{\infty}\mathrm{d}\tau p(\tau)$ where $p(t)$ is the probability
density for observing the event at time $t$ [2]. Perhaps the most widely
adopted model is Cox regression [5] which fits the so-called hazard function
$h(t)\equiv p(t)/S(t)$ via a simple (log-) linear function on top of a
baseline hazard $h_{0}(t)$ as $h(t)=h_{0}(t)\exp[\bm{\beta}\cdot\bm{\bm{x}}]$
with weights $\bm{\beta}$.

In addition, there are also efforts to combine machine learning techniques
with survival analysis. For instance, staying within Cox’s [5] proportional
hazards (PH) [i.e., $h(t)\propto\exp(\bm{\beta}\cdot\bm{\bm{x}})$] setting,
one can use boosting [6] to learn the parameters, or extend the linear
function using a neural network architecture, as done in Refs. [7, 8, 9].
Neural network structures that go beyond the PH assumption usually rely on the
binning of the time-to-event variables [10, 11, 12]. Apart from neural
networks, other models such as random forests [13] and support vector machines
[14, 15, 16, 17] have been developed as well.

While these models focus on one survival variable, joint analysis of multiple
unordered and distinct time recordings has received comparitavely little
attention. A more traditional statistical approach, such as the Wei-Lin-
Weissfeld model [18] (see also Ref. [19] for a related review), solves a Cox
model for each individual time-to-event variable and subsequently performs
joint inference on the parameters to determine their significance. MEPSUM [20]
is a more recent mixture model where each mixing component fits a discrete
time hazard function.

In this work, a different approach is taken by (unsupervised) training an
energy based model on the multiple time-to-event variable likelihood function
[21]. Specifically, we consider a neural network inspired model called a
harmonium [22] (or restricted Boltzmann machine [23], as it is synonymously
called), and adapt it to the approach of survival analysis.

This paper is organised as follows: first, the theory behind the model,
together with a training procedure, is laid out in the 2 section. Next, we
illustrate that a harmonium can capture non-linear survival patterns. We
subsequently benchmark the harmonium to established off-the-shelf survival
models in the 4. Finally, we discuss the relevance and scope of our findings
in the 5.

##  2 Theory

Before building a new harmonium for survival data, let us first briefly review
some of its quintessential properties.

###  2.1 Harmoniums

The textbook harmonium [24] consists of binary input states $\bm{x}$
($x_{i}\in\\{0,1\\}$ where $i=1\dots n_{v}$), binary activations $\bm{h}$
($h_{j}\in\\{0,1\\}$ where $j=1\dots n_{h}$), and an energy function $E$ that
linearly couples $\bm{x}$ and $\bm{h}$ through a receptive field $\bm{W}$ as

$E(\bm{x},\bm{h})=\bm{x}^{T}\bm{W}\bm{h}.$ (1)

The energy function encodes the preference for specific assignments of
$\bm{x}$ and $\bm{h}$ (with low $E$) over other assignments for which $E$ is
higher.

The probability distribution is parametrised by the energy $E$ as

$p(\bm{x},\bm{h})=\frac{1}{Z}e^{-E(\bm{x},\bm{h})},$ (2)

and a normalisation constant $Z$, called the partition function, which only
depends on the adjustable parameters (in this case $\bm{W}$). The latent
states themselves are not observed but enrich
$p(\bm{x})=\sum_{\bm{h}}p(\bm{x},\bm{h})$’s capacity to capture higher-order
(i.e., beyond pair-wise) statistics in the data [25]. Unfortunately, since the
partition function $Z$ is intractable [26] so is $p(\bm{x})$. Sampling from
$p(\bm{x}|\bm{h})$ and $p(\bm{h}|\bm{x})$ is nevertheless easy thanks to the
bipartite structure of $E(\bm{x},\bm{h})$. The interpretation of $\bm{W}$ as a
receptive field derives from the activation function of $h_{j}$ given the
visible states $\bm{x}$, i.e.,
$p(h_{j}=1|\bm{x})=\sigma(-\sum_{i=1}^{n_{v}}x_{i}W_{ij})$ with sigmoid
activation function $\sigma(x)=1/(1+\exp[-x])$, which is structurally akin to
a neural network. Likewise, the latent states activate the visible binary
states $\bm{x}$ through
$p(x_{i}=1|\bm{h})=\sigma(-\sum_{j=1}^{n_{h}}W_{ij}h_{j})$.

Given a set of $m$ samples $\\{\bm{x}^{(i)}\\}_{i=1}^{m}$, training proceeds
by adjusting the free parameters $\bm{\Theta}$ in $E(\bm{x},\bm{h})$—which in
this case consists of $\bm{W}$—to maximise the log-likelihood function

$\mathcal{L}(\\{\bm{x}^{(i)}\\}_{i=1}^{m})=\frac{1}{m}\sum_{i=1}^{m}\ln
p(\bm{x}^{(i)}),$ (3)

by approximating its gradient using Gibbs samples. The contrastive divergence
algorithm relies on the decomposition of the free parameter $\bm{\Theta}$
gradient of the likelihood

$\nabla_{\bm{\Theta}}\mathcal{L}=-\left(\left\langle\nabla_{\bm{\Theta}}E\right\rangle_{p(\bm{h}|\bm{x})p_{\mathrm{data}}(\bm{x})}-\left\langle\nabla_{\bm{\Theta}}E\right\rangle_{p(\bm{x},\bm{h})}\right),$
(4)

into an expectation over the empirical data [first term on the right hand side
(rhs),
$p_{\mathrm{data}}(\bm{x})=\frac{1}{m}\sum_{i=1}^{m}\delta_{\bm{x},\bm{x}^{(i)}}$]
called the positive phase and an expectation over the model itself
($\nabla_{\bm{\Theta}}\ln
Z=-\langle\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\rangle_{p(\bm{x},\bm{h})}$,
second term rhs) referred to as the negative phase [24].

While the positive phase can be evaluated in closed form, the negative phase
(i.e., the partition function gradient [24]
$-\langle\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\rangle_{p(\bm{x},\bm{h})}$) is
to be approximated by Gibbs sampling between $p(\bm{x}|\bm{h})$ and
$p(\bm{h}|\bm{x})$. The key emperical observation behind the contrastive
divergence algorithm [23, 27] is that, when initialising the chain with the
data distribution, a single chain step usually suffices to estimate the
negative phase.

###  2.2 Model

Having reviewed the basic structure of harmoniums, lets now turn to the
survival data. To accommodate survival data we need to: (i) design an energy
function capturing both survival, categorical and other numerically valued
variables (Sec. 2.2.1), (ii) adapt the objective (i.e., likelihood) function
to account for censoring and completely missing data (Sec. 2.2.2), and (iii)
layout a corresponding training algorithm (Sec. 2.2.3). Let us first focus on
the energy function.

####  2.2.1 Energy function

Figure 1: Graphical representation of the energy function $E(\bm{x},\bm{h})$.
Four types of variables (nodes) can be distinguished: binary states
$\bm{x}^{A}$, time-to-event variables $\bm{x}^{B}$, continuous variables
$\bm{x}^{C}$, and (unobserved) binary latent states $\bm{h}$. Edges between
$\bm{x}$ and $\bm{h}$, indicating the receptive fields, form a bipartite
graph.

To reiterate, the energy function $E$ codifies preferences for specific
assignments of the variables. Since the event times are numerically valued
(instead of binary, like in Sec. 2.1), we adjust the energy function
accordingly. Beside the the survival events, we set out to capture binary
variables (e.g., smoking status) and numeric factors (like body mass index).
Let us therefore distinguish between three sets of input variables (denoted by
$A$, $B$, and $C$):

  * 1.

Categorical variables $\bm{x}^{A}=\\{x_{i}:i\in A\\}$ that are encoded as
binary variables $x_{i\in A}\in\\{0,1\\}$.

  * 2.

Time-to-event variables $\bm{x}^{B}=\\{x_{i}:i\in B\\}$ that are appropriately
scaled to the unit interval $(0,1]$.

  * 3.

Other numerical factors $\bm{x}^{C}=\\{x_{i}:i\in C\\}$ that can take on
arbitrary values on the real line $x_{i\in C}\in\mathcal{R}$.

Both the binary variables (in $A$) and the numeric factors (in $C$) are
assumed to be independent of time. (The interval on which the time-to-event
variables are defined is discussed in more detail in A.) In addition, a single
set of latent variables $\\{h_{i}\\}_{i\in H}$ is used (like in Sec. 2.1) to
coherently model the three sets of input variables. The latent code will be
restricted to binary values $h_{i}\in\\{0,1\\}$ in view of its regularising
effect [27].

Assigning an energy term to each variable type gives rise to the overall
energy function

$E(\bm{x},\bm{h})=E_{A}(\bm{x}^{A},\bm{h})+E_{B}(\bm{x}^{B},\bm{h})+E_{C}(\bm{x}^{C},\bm{h})+E_{H}(\bm{h}).$
(5)

Observe that the energy function’s bi-partite structure between $\bm{x}$ and
$\bm{h}$ determines the form of $p(\bm{x}|\bm{h})$. Decomposing $p(\bm{x})$ as
$\sum_{\bm{h}}p(\bm{x}|\bm{h})p(\bm{h})$ reveals that the conditional
distributions $p(\bm{x}|\bm{h})$ are the building blocks of the model,
weighted by the ($\bm{\Theta}$ dependent) latent priors $p(\bm{h})$. So the
choice of the energy term controls the building blocks of the model.

For the categoric $\bm{x}^{A}$ and numeric factors $\bm{x}^{C}$ we rely on
established energy functions: (i) $E_{A}$ is modelled as a binary-binary
harmonium [27], similar to the energy function Sec. 2.1 but supplemented with
a bias term, (ii) $E_{C}$ represents a Gaussian-binary harmonium [28, 27, 29].

For the time-to-event variables a new function $E_{B}$ is proposed (detailed
in A) for which $p(x_{i}|\bm{h})$ is a (truncated) gamma distribution and use
$E_{H}$ to bias the latent states $\bm{h}$. In summary, we take the following
forms for the individual energy functions

$\displaystyle E_{A}$ $\displaystyle=$ $\displaystyle\sum_{i\in A,j\in
H}x_{i}W^{A}_{ij}h_{j}+\sum_{i\in A}x_{i}a^{A}_{i}\,,$ (6) $\displaystyle
E_{B}$ $\displaystyle=$ $\displaystyle\sum_{i\in B,j\in
H}x_{i}W^{B}_{ij}h_{j}-\ln(x_{i})\|V_{ij}\|h_{j}+\sum_{i\in
B}x_{i}a^{B}_{i}-\ln(x_{i})\|c_{i}\|\,,$ (7) $\displaystyle E_{C}$
$\displaystyle=$ $\displaystyle\sum_{i\in C,j\in
H}\frac{x_{i}W^{C}_{ij}h_{j}}{\sigma_{i}}+\sum_{i\in
C}\frac{(x_{i}-a^{C}_{i})^{2}}{2\sigma_{i}^{2}}\,,$ (8) $\displaystyle E_{H}$
$\displaystyle=$ $\displaystyle\sum_{j\in H}b_{j}h_{j}\,.$ (9)

The weights $\\{\bm{W}^{A},\bm{W}^{B},\bm{W}^{C}\\}$ can be interpreted as the
receptive fields of $\bm{x}$ to activate the latent states $\bm{h}$, while
$\\{\bm{a}^{A},\bm{a}^{B},\bm{a}^{C}\\}$ and $\bm{b}$ are their respective
biases. The $\bm{V}$ and $\bm{c}$ terms in $E_{B}$ are additional receptive
fields and biases that help modulate the survival distribution. The receptive
fields of $E(\bm{x},\bm{h})$, coupling $\bm{x}$ and $\bm{h}$, are illustrated
in Fig. 1 by corresponding edges. The form of $E(\bm{x},\bm{h})$ fixes the
distribution over $\bm{x}$ given $\bm{h}$, and it is not difficult to show
that

$p(x_{i}|\bm{h})=\left\\{\begin{matrix}\sigma[(1-2x_{i})z_{i}]&i\in A,\\\
p_{\Gamma}(x_{i}|\alpha_{i},\beta_{i})&i\in B,\\\
\mathcal{N}(x_{i}|\mu_{i},\sigma_{i}^{2})&i\in C.\end{matrix}\right.$ (10)

where $\sigma[(1-2x_{i})z_{i}]$ is the sigmoid function with latent state
activation $z_{i}=a^{A}_{i}+\sum_{j\in H}W^{A}_{ij}h_{j}$,
$p_{\Gamma}(x_{i}|\alpha_{i},\beta_{i})$ the right truncated Gamma
distribution [Eq. (21)] with shape $\alpha_{i}=\sum_{j\in
H}\|V_{ij}\|h_{j}+\|c_{i}\|+1$ and rate $\beta_{i}=\sum_{j\in
H}W^{B}_{ij}h_{j}+a^{B}_{i}$, and finally
$\mathcal{N}(x_{i}|\mu_{i},\sigma_{i}^{2})$ a Gaussian with mean
$\mu_{i}=a^{C}_{i}-\sigma_{i}\sum_{j\in H}W^{C}_{ij}h_{j}$ and standard
deviation $\sigma_{i}$. In a similar way, the activations of the latent
variables

$p(h_{j}|\bm{x})=\sigma[(1-2h_{j})\phi_{j}],$ (11)

depend on the contributions of all variable types, which are jointly captured
by

$\phi_{j}=b_{j}+\sum_{i\in A}x_{i}W^{A}_{ij}+\sum_{i\in
B}x_{i}W^{B}_{ij}+\sum_{i\in C}\frac{x_{i}W^{C}_{ij}}{\sigma_{i}}-\sum_{i\in
B}\ln(x_{i})\|V_{ij}\|.$ (12)

As explained in Sec. 2.1, a key observation that is central to the training of
harmoniums is that Gibbs samples from $p(\bm{x},\bm{h})$ can be obtained by
alternating between Eq. (10) and Eq. (11). In this way, an entire block of
states can be updated in parallel, thanks to its conditional independence.

####  2.2.2 Cost objective

With the energy function in place, let us now adapt the likelihood function
(our cost objective) to incorporate partially and completely missing values.
We will assume uninformative censoring and that missing values are missing at
random. To simplify the exposition we shall henceforth focus on right censored
data (or censored, for short). That is, observations for which there is a
lower bound on the failure time (e.g., a participant that was lost to follow-
up after time $t$). The standard likelihood approach for modelling a single
time-to-event variable is as follows. In the presence of censoring, the
likelihood is modified by replacing $p(t)$ by its survival function $S(t)$
when it is censored at time $t$. Writing
$S(t)=\int_{0}^{\infty}\mathrm{d}\tau\Theta(\tau-t)p(\tau)$ with the Heaviside
step function $\Theta(x)=\left\\{\begin{matrix}1&x\geq 0\\\
0&x<0\end{matrix}\right.,$ reveals that, in fact, we are integrating out the
unobserved region. Both cases, censored and observed events, can be succinctly
codified as an integration over the domain of event times
$\int_{0}^{\infty}\mathrm{d}\tau p(\tau)\chi(t,\tau,e)$ constrained by
$\chi(t,\tau,e)=\delta(\tau-t)^{e}\Theta(\tau-t)^{1-e}$ with $\delta(x)$ the
Dirac delta function and $e=1$ ($e=0$) indicating observation (censoring) at
time $t$. Note the analogy with completely missing data where the entire
domain (instead of subset of the domain) is marginalised, e.g.,
$p(t_{1})=\int_{-\infty}^{\infty}\mathrm{d}t_{2}p(t_{1},t_{2})$ when
$t_{2}\in\mathcal{R}$ is missing. We can therefore apply a similar
codification scheme to missing values to obtain
$\chi(t,\tau,e)=\delta(\tau-t)^{e}$. A consistent generalisation from one to
multiple censored variables is straightforward: integrate out the entire
unobserved region [21]. To make this more precise, let $e_{a}=1$ indicate the
occurrence and $e_{a}=0$ the absence of observation $x_{a}$. That is,
$e_{a}=0$ indicates that $x_{a}$ is censored when $a$ refers to a time-to-
event variable ($a\in B$), or completely missing otherwise ($a\in A\cup C$).
In addition, let $\xi_{a}$ be the corresponding observed value when $e_{a}=1$,
its lower bound (i.e., censoring time) for the survival variables (i.e., $a\in
B$) or a placeholder $\xi_{a}=?$ otherwise ($a\in A\cup C$) when $e_{a}=0$. As
a shorthand, denote $o_{a}=(\xi_{a},e_{a})$ and a superscript $o_{a}^{(i)}$ to
refer to a specific sample $i$. First, group the marginalisation constraints
that are imposed by the observations

$\chi(\bm{x},\bm{o})=\prod_{i\in A\cup B\cup
C}\delta(x_{i}-\xi_{i})^{e_{i}}\prod_{j\in B}\Theta(x_{j}-\xi_{j})^{1-e_{j}},$
(13)

with $\Theta(x)$ the Heaviside step function and $\delta(x)$ the Dirac delta
function (Kronecker delta function) for the continuous variables in $B\cup C$
(binary variables in $A$). Equation (13) is a symbolic way to represent that
we should either pick the observed values (the delta function) or marginalise
the unobserved region (which, for the survival variables, is the interval
starting from the censor time, or the entire domain otherwise). In this way,
the likelihood can be expressed as

$L(\\{\bm{o}^{(i)}\\}_{i=1}^{m})=\prod_{i=1}^{m}\int\mathrm{d}\bm{x}p(\bm{x})\chi(\bm{x},\bm{o}^{(i)}),$
(14)

using the shorthand
$\int\mathrm{d}\bm{x}\equiv\int_{-\infty}^{\infty}\mathrm{d}\bm{x}^{C}\int_{0}^{1}\mathrm{d}\bm{x}^{B}\sum_{\bm{x}^{A}\in\\{0,1\\}\otimes^{|A|}}$.
In summary, to train the model that takes into account censored and missing
data should strive to optimise the likelihood function Eq. (14) or,
equivalently, the log-likelihood function
$\mathcal{L}(\\{\bm{o}^{(i)}\\}_{i=1}^{m})=\ln(L)/m$.

Having spelled out the likelihood function in fair generality, next we apply
it to the energy parameterisation
$p(\bm{x})\propto\sum_{\bm{h}}\exp[-E(\bm{x},\bm{h})]$. To keep the bi-partite
structure of the model for sampling [which would be destroyed by directly
carrying out the integration in Eq. (14)] we turn to a trick from Ref. [30] to
reformulate the model in terms of $p(\bm{o},\bm{x},\bm{h})$ [30], where

$p(\bm{o},\bm{x},\bm{h})\propto e^{-E(\bm{x},\bm{h})}\chi(\bm{x},\bm{o}).$
(15)

With the help of Eq. (15) the gradient of the log-likelihood (details are in
B) can be expressed as

$\nabla_{\bm{\Theta}}\mathcal{L}=-\left(\frac{1}{m}\sum_{i=1}^{m}\left\langle\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\right\rangle_{p(\bm{x},\bm{h}|\bm{o}^{(i)})}-\langle\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\rangle_{p(\bm{x},\bm{h})}\right).$
(16)

Heuristically speaking, Eq. (16) indicates that the gradient contrasts the the
empirical statistics of $\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})$ [through
$p(\bm{x},\bm{h}|\bm{o}^{(i)})$, first term, rhs] with the models own
perception [generated by $p(\bm{x},\bm{h})$] of
$\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})$ (second term, rhs).

####  2.2.3 Training

1:while not converged do

2: Load minibatch $\\{\bm{o}^{(1)},\dots,\bm{o}^{(m)}\\}$ of $m$ samples.

3: for $i=1$ to $m$ do

4: $\bm{x}^{(i)}\leftarrow\bm{\xi}^{(i)}$

5: $\tilde{\bm{x}}^{(i)}\leftarrow\bm{\xi}^{(i)}$

6: for $l=1$ to $k$ do

7: $\bm{h}^{(i)}\sim p(\bm{h}^{(i)}|\bm{x}^{(i)})$ $\triangleright$ Sample
positive phase.

8: $\bm{x}^{(i)}\sim p(\bm{x}^{(i)}|\bm{h}^{(i)},\bm{o}^{(i)})$

9: $\tilde{\bm{h}}^{(i)}\sim p(\tilde{\bm{h}}^{(i)}|\tilde{\bm{x}}^{(i)})$
$\triangleright$ Sample negative phase.

10: $\tilde{\bm{x}}^{(i)}\sim p(\tilde{\bm{x}}^{(i)}|\tilde{\bm{h}}^{(i)})$

11: end for

12: $\bm{\mu}^{(i)}\leftarrow p(\bm{h}^{(i)}=\bm{1}|\bm{x}^{(i)})$

13: $\tilde{\bm{\mu}}^{(i)}\leftarrow
p(\tilde{\bm{h}}^{(i)}=\bm{1}|\tilde{\bm{x}}^{(i)})$

14: end for

15: $\triangleright$ Gradient ascent update with learning rate
$r_{\mathrm{learn}}$.

16:
$\bm{\Theta}\leftarrow\bm{\Theta}-\frac{r_{\mathrm{learn}}}{m}\sum_{i=1}^{m}\left(\nabla_{\bm{\Theta}}E(\bm{x}^{(i)},\bm{\mu}^{(i)})-\nabla_{\bm{\Theta}}E(\tilde{\bm{x}}^{(i)},\tilde{\bm{\mu}}^{(i)})\right)$.

17:end while

Algorithm 1 The $k$-step contrastive divergence algorithm for censored and
missing values.

Next, we discuss how to approximate the gradient [Eq. (16)] in order to
maximise the likelihood with gradient ascent. In the standard contrastive
divergence [23] approach, the negative phase is approximated using Gibbs
samples while the positive phase can be evaluated exactly. Incorporating
missing and censored values has modified the positive phase [first term, rhs
Eq. (16)] in a way that evades a closed form solution. Instead, Eq. (16) will
be estimated by Gibbs sampling both phases.

Analogous to the sampling of $p(\bm{x},\bm{h})$ for the negative phase, we
alternate between $p(\bm{x}|\bm{o},\bm{h})$ and $p(\bm{h}|\bm{o},\bm{x})$ to
generate $p(\bm{x},\bm{h}|\bm{o})$ samples for the positive phase where

$p(x_{i}|\bm{o},\bm{h})=\left\\{\begin{matrix}\delta(x_{i}-\xi_{i})&\forall
i&e_{i}=1,\\\ p(x_{i}|\bm{h})&i\in A\cup C&e_{i}=0,\\\
p^{\Gamma}_{[\xi_{i},1]}\left[x_{i}\left|\alpha_{i}(\bm{h}),\beta_{i}(\bm{h})\right.\right]&i\in
B&e_{i}=0,\\\ \end{matrix}\right.$ (17)

with
$p^{\Gamma}_{[\xi_{i},1]}\left[x_{i}|\alpha_{i}(\bm{h}),\beta_{i}(\bm{h})\right]$
the gamma distribution normalised to the $[\xi_{i},1]$ interval [Eq. (43)] and

$p(\bm{h}|\bm{o},\bm{x})=p(\bm{h}|\bm{x}).$ (18)

Physically, Eq. (17) indicates that samples $x_{i}$ should adhere to the
bounds imposed by the observation (which is the lower-bound censor time for
the censored variables).

In summary, the training procedure is as follows: For the positive phase
[first term r.h.s. Eq. (16)] replace the censored and missing values with
“fantasy states” and calculate the $\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})$
statistic. To replace the censored data, clamp the observed events and sample
the censored events from the unobserved interval and sample the missing values
from the entire distribution [Eq. (17)]. The negative phase [second term
r.h.s. Eq. (16)] is calculated similarly, but now all the states are updated
(none are clamped) from the entire interval [Eq. (10)]. To generate fantasy
states, pick the the mini batch as the initial state of the Gibbs chain,
initialise the placeholder “$?$” by the median value over the training set,
and carry out $k$ Gibbs chain steps. Finally, update the weights using the
phase difference and repeat the entire process until some predefined stopping
criterion. In pseudocode, the algorithm is outlined in Algorithm 1.
Reassuringly, the original contrastive divergence algorithm is recovered as a
special case when all the data is observed (i.e., $e^{(i)}_{a}=1$ for all $i$
and $a$).

##  3 Example: a non-linear two-dimensional survival distribution

(a)

(b)

Figure 2: A harmonium captures non-linear related time recordings. (a) A
synthetic time-to-event distribution where the timing of two separate events,
$t_{1}$ and $t_{2}$, are non-linearly related through colour. The density is
composed of two red ($v^{(1)}$ and $v^{(2)}$, dashed contours) and two blue
($v^{(3)}$ and $v^{(4)}$, solid contours) gamma distributed blobs. Side panels
show the marginal probability density by colour, illustrating the multivariate
nature of the problem. (b) Model fit (contours) of observed (dots) and
censored recordings (crosses) sampled from (a). The contours, indicating
constant probability density of the harmonium, shows that all four colour-mode
combinations are recapitulated.

The limitations of uni-survival variate (i.e., conventional survival) models
is best illustrated with an intrinsically two-dimensional problem. Consider a
data generating distribution yielding two event recordings
$\bm{x}^{B}=[t_{1},t_{2}]^{T}$ and a colour, red ($x^{A}=0$, we’ve dropped the
index for convenience) or blue ($x^{A}=1$). Let the probability density be
confined to the unit square $[0,1]\times[0,1]$ symmetrically tiled with four
equally weighted bell-shaped blobs: blue generating recordings along the
diagonal, and red on the off-diagonal quadrants (see Fig. 2(a)). Specifically,
blobs are two-dimensional, independent (i.e., $t_{1}\perp t_{2}$), unit-
interval truncated Gamma distributions [Eq. (21)] with modes placed at
$\bm{v}^{(1)}=\left(\frac{1}{4},\frac{1}{4}\right)$,
$\bm{v}^{(2)}=\left(\frac{3}{4},\frac{3}{4}\right)$ corresponding to red
($x^{A}=0$) and $\bm{v}^{(3)}=\left(\frac{3}{4},\frac{1}{4}\right)$,
$\bm{v}^{(4)}=\left(\frac{1}{4},\frac{3}{4}\right)$ for blue ($x^{A}=1$).
Modes are sufficiently squeezed (shape and rate $\alpha=8.1$, $\beta=58$ or
$\alpha=29$, $\beta=76$) so as to form a Gaussian-like shape and sampled with
equal probability.

Looking at the projections (i.e., marginals) along the axes (side panels Fig.
2(a)), shows how the red and blue modes collapse onto each other. Viewed from
either $t_{1}$ or $t_{2}$ alone, one would therefore be inclined to (falsely)
conclude there is no relation between colour and survival.

On the other hand, we can derive an approximate, not necessarily unique, but
nevertheless closed-form, solution for a harmonium with $n_{h}=4$ hidden
states. The harmonium approximates a mixture of Gaussians in survival space
(parametrised by $\bm{x}^{B}$) with $x^{A}$ clamped to mode $j$’s colour
${\tilde{x}^{A}_{j}}$ (see C for a derivation). The probability density can be
approximated as

$p(x^{A},\bm{x}^{B})=\frac{1}{Z}\sum_{h_{1},\dots,h_{4}}\exp[-E(x^{A},\bm{x}^{B},\bm{h})]\approx\frac{1}{4}\sum_{j=1}^{4}\delta_{x^{A},\tilde{x}^{A}_{j}}\mathcal{N}\left(\bm{x}^{B}|\bm{v}^{(j)},\bm{\Sigma}_{j}\right),$
(19)

where the mean $\bm{v}^{(j)}$ and the (diagonal) covariance matrix
$\bm{\Sigma}_{j}$ are determined through the rows of the receptive fields
$\bm{v}^{(j)}=\frac{\|\bm{V}_{j}\|}{\bm{W}^{B}_{j}}$ and
$\bm{\Sigma}_{j}=\frac{\|\bm{V}_{j}\|}{(\bm{W}^{B}_{j})^{2}}$, with all other
visible biases zero (all other weights are described C).

Knowing that the problem is (approximately) solvable in theory, lets
illustrate training when some of the time recordings are censored. We used the
aforementioned distribution to generate 1000 samples and censored each event
time $t_{i}>\frac{3}{4}$ with 75 % probability. For clarity, half of the
points are shown in Fig. 2(b), coloured by $x^{A}$, and marked by a cross
where censored. The harmonium was trained for $3\cdot 10^{5}$ epochs with a
learning rate $r_{\mathrm{learn}}=0.375$, 10 % momentum, and 3 persistent [31]
contrastive divergence sampling steps.

While a model that doesn’t account for censoring would underestimate survival,
we find that the harmonium correctly identifies all four modes. We do observe
that the modes are less sharply peaked (more smeared) compared to the original
distribution. This is attributed to the approximate and stochastic nature of
the contrastive divergence algorithm, which sometimes hinders convergence.
Overall, the harmonium satisfactory captures the non-linear relationship in
the survival data. As a reference, we trained a Cox model [5] on either
$t_{1}$ or $t_{2}$ with $x^{A}$ as a covariate. In both cases we found that
its regression coefficient is zero (null hypothesis) under a p-value threshold
of 0.05. That is, the Cox model finds no relation between $x^{A}$ and
survival.

##  4 Experiments

To illustrate performance on real world datasets, the harmonium is compared to
(i) Cox regression [5] from the lifelines package [32] with both $L_{1}$ and
$L_{2}$ regularisation, (ii) random survival forest [13] and (iii) the fast
support vector machine (SVM) [17], where the latter two are both from the
scikit-survival package.

###  4.1 Datasets

Our benchmark is comprised of four lifelines datasets [32], namely:

  * 1.

The recidivism of convicts released from the Maryland state prisons ($m$=432
convicts) [33]—denoted as arrest—to study the effect of financial aid. All
features were used for training, but we grouped the number of prior
convictions $>$ 5 before one-hot encoding categories as dummies.

  * 2.

The duration of democratic and dictatorial political regimes ($m$=1808
countries) [34]—denoted as democracy—with the continent name and type of
regimes as features.

  * 3.

The survival of women with breast cancer ($m$=686 patients) [35] (denoted as
gbsg2) to measure the effect of hormonal therapy, using all features.

  * 4.

The survival of advanced lung cancer patients ($m$=288 patients) [36]—denoted
as ncctg—where the prognostic value of a patient’s questionnaire was examined.
We used all features except for the institute code and the columns weight loss
and meal calory intake. After one-hot encoding caterories as dummies, the low
variance features that were on/off in more than 95 % of the samples were
dropped to prevent collinearity.

In addition, our benchmark comprises two additional lung cancer datasets
containing two time-to-event variables (bundled with the code, but not part of
lifelines):

  * 1.

A Dutch study, nvalt11, considered the effect of profylactic brain radiation
versus observation in ($m$=174) patients with advanced non-small cell lung
cancer [37]. The nvalt11 dataset contained time recordings for overall
survival (OS) and symptomatic brain metastasis-free survival (SBMFS). We
modelled the categoric variables gender, control arm, performance status, and
smoking status plus the numeric feature age. The categorical features
histology, prior medical conditions, prior malignancies, and stage and numeric
variable BMI contained missing values, and were therefore dropped in all
models except the harmonium.

  * 2.

Another Dutch study, called nvalt8 ($m$=200 participants), that examined if
nadroparin combined with chemotherapy could reduce cancer relapse after
surgical removal of a non-small cell lung tumour [38]. The dataset contained
failure times for both OS and recurrence free survival (RFS). All models used
the numeric feature age and the categories: performance status, histology,
smoking status, stage of the disease, control arm and the T, N, and M tumour
classification categories. After one-hot-encoding, the low variance features
that were on/off in more than 95 % of the samples were dropped. The harmonium
also incorporated the metabolic activity measured as FDG-PET SUVmax $\geq$10
and the numeric variable BMI that both contained missing values.

Apart from the time-to-event variables, all numeric features were standardised
and categorical variables were dummy encoded prior to training.

###  4.2 Results

(a) Discriminative accuracy (higher is better).

(b) Calibration (lower is better).

Figure 3: A benchmark of survival models (indicated in the legend) across
various datasets (horizontal axis) shows that the harmonium is on par with
uni-survival variate models, and discriminative predictions improve with an
additional survival variable (nvalt11 and nvalt8). Specically, overall
survival (OS) and symptomatic brain metastasis-free survival (SBMFS),
$\bm{x}^{B}=[x_{\mathrm{OS}},x_{\mathrm{SBMFS}}]^{T}$, were recorded for the
nvalt11 dataset; OS and recurrence free survival (RFS),
$\bm{x}^{B}=[x_{\mathrm{OS}},x_{\mathrm{RFS}}]^{T}$, were recorded for the
nvalt8 dataset. Metrics for these two datasets were computed for OS but we
distinguish between two methods of computation for the harmonium. Namely,
metrics that factor in the second survival variable [indicated by a * in the
legend, and computed through survival distribution
$S(x_{\mathrm{OS}}|\bm{o}_{-\mathrm{OS}})$] versus the metrics where this
variable was marginalised out (without a *, via
$S(x_{\mathrm{OS}}|\bm{o}_{-\\{OS,SBMFS/RFS\\}})$). Calibration data was not
available for the support vector machine (svm). Markers and errorbars indicate
the mean value and the standard deviation from 5x5 nested cross-validation.

To reiterate, the primary difference between the harmonium, and the
implementations of Cox model, random survival forest, and SVM considered here,
is that the harmonium can incorporate missing values and multiple (potentially
non-linearly related) survival variables. Results are compared using two
metrics: Harrell’s concordance index [1] and Brier’s calibration loss [39]
(details are in D.1).

Since the harmonium captures both survival recordings of the nvalt datasets
and computes the metrics based on the survival distribution (see Sec.D.1), we
can compute it in two ways. Firstly, by factoring in the additional survival
variable $S(x_{\mathrm{OS}}|\bm{o}_{-\mathrm{OS}})$ where
$\bm{o}_{-\mathrm{OS}}$ denotes observation $\bm{o}$ with the element
$\mathrm{OS}$ removed (but still containing the other survival variable).
Alternatively, the dependence of the other survival variable (SBMFS and RFS)
can be marginalised out giving
$S(x_{\mathrm{OS}}|\bm{o}_{-\\{\mathrm{OS},\mathrm{SBMFS}\\}})$ and
$S(x_{\mathrm{OS}}|\bm{o}_{-\\{\mathrm{OS},\mathrm{RFS}\\}})$ for nvalt11 and
nvalt8, respectively.

The benchmark results, summarised in Fig. 3, shows comparable performance
across all models for the first four datasets. For the nvalt datasets, notice
that when we factor in the additional survival information (indicated by a *
in the legend) we observe a substantial improvement in the model’s
predictions, especially in terms of the concordance index. For the nvalt8
dataset, we find that the harmonium does slightly worse in terms of the
calibration (Fig. 3(b)), but overall scores comparable over the range of
datasets. These results are in line with common sense, that a disease relapse
or finding a brain tumour decreases one’s expected life expectancy.
Conversely, when marginalising out the extra time-to-event variable for the
nvalt datasets (without a *) the performance reduces to that of the other
models. This performance reduction upon marginalisation is illustrative for
the information content captured by the complementary endpoint. Including
variables with missing values, as we did for the nvalt datasets, showed no
noticeable improvement for the harmonium compared to the other models (where
this could not be taken into account).

##  5 Discussion & Conclusion

Healthcare data follows an inherent timeline where new information, such as a
lab result or a diagnosis, comes in continuously. At the same time, both
statistical and machine learning models require data in tabular format. This
poses a challenge, where one should strike a balance between a format that
accommodates the model and simultaneously does justice to the time ordering of
the data. Our work is a step towards a consolidation of these two
representations of the data, by modelling both missing values and multiple
time-to-event variables in one coherent framework.

One important disadvantage of the present unsupervised model is that, like
many neural networks, there is a myriad of hyperparameters to tune. Choosing
appropriate parameters for the learning rate, batch size, number of latent
states, how many epochs to train, and regularisation can be challenging. In
addition, while some quantities, e.g., the latent states, can be computed
efficiently (i.e., linear in the number of input variables), others such as
the survival distribution [Eq. (42) in D] can be computationally more
demanding.

A second limitation of this work, unrelated to our model, is that Harrell’s
concordance index and Brier loss — both intrinsically uni-survival variate
metrics — may not be the most appropriate measures to comprehensively
interrogate a model’s capacity to capture multiple time recordings. We could
only indirectly probe its performance by conditioning on, and marginalising
out, the second survival variable. Alas, as far as we know, no higher
dimensional generalisations of, e.g., Harrell’s concordance index or the Brier
loss exist. In this regard, we believe that the two-dimensional survival
distribution in Sec. 3 can serve as a useful example, helping future
researchers test their models on multi-dimensional non-linear survival data.

In summary, a new harmonium was proposed for partially and completely missing
data. Multiple distinct failure times can be jointly modelled without imposing
a priori ordering or relation amongst the events, in contrast to conventional
survival techniques. In addition, time-independent categorical and numerical
features with missing values can be straightforwardly incorporated thanks to
its generative structure. Our work shows that the harmonium can capture non-
linear survival patterns and illustrates the information embodied in
complementary survival endpoints while accounting for the uncertainty in the
failure time. We have eliminated the need for selecting a single endpoint and
pave the way towards a unified timeline view of the data.

## Acknowledgements

We thank Rik Huijzer for proofreading the manuscript. We would like to thank
the Center for Information Technology of the University of Groningen for
providing access to the Peregrine high performance computing cluster.

## Competing interests

The authors declare that they have no competing interests.

## Availability of data and materials

Code, data, and examples are publicly available under the open source Apache
License 2.0 at https://gitlab.com/hylkedonker/harmonium-models.

##  Appendix A Energy function survival variables

All survival variables are assumed to be offset against a fixed landmark
(e.g., the date when a participant entered the study) so that all values are
$>0$ and we focus on right censored events. Our goal is to construct an energy
function for the survival variables where the building blocks are composed of
gamma distributions. To this end, we make the following ansatz for the energy
function

$E_{B}(\bm{x}^{B},\bm{h})=\sum_{i\in B}\sum_{j\in
H}x_{i}W^{B}_{ij}h_{j}-\ln(x_{i})\|V_{ij}\|h_{j}+\sum_{i\in
B}x_{i}a^{B}_{i}-\ln(x_{i})\|c_{i}\|.$ (20)

Henceforth $x_{i}$ are assumed to be scaled to the unit interval
$x_{i}\in(0,1]$ for all $i\in B$ by normalising $x_{i}\rightarrow
x_{i}/\tau_{i}$ with a suitably chosen time horizon $\tau_{i}$. Our motivation
is two fold. The first is technical: $x_{i}\in(0,1]$ ensures that the
probability density functions [Eq. (21), below] can be normalised without
imposing the constraint $\beta_{i}>0$ for $\beta_{i}=\sum_{j\in
H}W^{B}_{ij}h_{j}+a^{B}_{i}$ (to prevent diverging integrals near infinity).
The second is physical: normalising the data reflects one’s believe about the
possible values the data can take since censoring at $\xi_{i}$ means: the
actual event is supposed to occur somewhere in the interval $[\xi_{i},1]$. For
example, it may be unrealistic to assume that a person becomes over 120 years
old, and the time horizon $\tau_{i}$ is a way to factor in these physical
constraints.

Invoking the definition of the model
$p(\bm{x}^{B},\bm{h})\propto\exp[-E_{B}(\bm{x}^{B},\bm{h})]$ yields the right
truncated gamma distribution

$p_{\Gamma}\left(x_{i}|\bm{h}\right)=\frac{x_{i}^{\alpha_{i}-1}e^{-\beta_{i}x_{i}}}{\Gamma(\alpha_{i})\gamma^{*}(\alpha_{i},\beta_{i})},$
(21)

as the conditional probability distribution, where $\alpha_{i}=\sum_{j\in
H}\|V_{ij}\|h_{j}+\|c_{i}\|+1$, $\Gamma(x)$ the Gamma function, and

$\gamma^{*}(a,z)=\frac{1}{\Gamma(a)}\int_{0}^{1}t^{a-1}e^{-zt}\mathrm{d}t,$
(22)

the incomplete gamma function [40]. Note that we haven’t exhausted the entire
parameter space by choosing the coupling strength $\|V_{ij}\|$ and bias
$\|c_{i}\|$ to be positive, whence $\alpha_{i}\geq 1$. Technically,
$\alpha_{i}$ must be larger than 0 to prevent poles from emerging. By
introducing a term $\|c_{i}\|\rightarrow\|c_{i}\|+d_{i}$ with $d_{i}>-1$
(e.g., $d_{i}=\lim_{A\uparrow 1}A\cos\varphi_{i}$) we can cover the entire
domain of $\alpha_{i}$. But to simplify the generation of samples from Eq.
(21) we focus on the form laid out in Eq. (20). Observe that the bias $\bm{c}$
can in principle be captured by $\bm{V}$ at the expense of introducing
additional latent states that are always turned on $h_{i}=1$. To reduce the
amount of parameters as much as possible, we choose instead to model the bias
$\bm{c}$ separate from $\bm{V}$.

##  Appendix B Derivation log-likelihood gradient

The derivation presented here parallels Ref. [30] with the appropriate changes
to the notation, and is provided here for completeness. Our goal is to
calculate the gradient of the log likelihood

$\mathcal{L}(\\{\bm{o}^{(i)}\\}_{i=1}^{m})=\frac{1}{m}\sum_{i=1}^{m}\ln\int\mathrm{d}\bm{x}p(\bm{x})\chi(\bm{x},\bm{o}^{(i)}),$
(23)

(abbreviated using
$\int\mathrm{d}\bm{x}\equiv\int_{-\infty}^{\infty}\mathrm{d}\bm{x}^{C}\int_{0}^{1}\mathrm{d}\bm{x}^{B}\sum_{\bm{x}^{A}\in\\{0,1\\}\otimes^{|A|}}$)
w.r.t. the expanded model

$p(\bm{o},\bm{x},\bm{h})=\frac{\exp[-E(\bm{x},\bm{h})]\chi(\bm{x},\bm{o})}{\Xi},$
(24)

where $\Xi$ is an unimportant normalisation constant which differs from the
partition function
$Z=\sum_{\bm{h}}\int\mathrm{d}\bm{x}\exp[-E(\bm{x},\bm{h})]$. To this end,
write
$Z(\bm{o})\equiv\sum_{\bm{h}}\int\mathrm{d}\bm{x}\exp[-E(\bm{x},\bm{h})]\chi(\bm{x},\bm{o})$
to further simplify the log likelihood to

$\mathcal{L}(\\{\bm{o}^{(i)}\\}_{i=1}^{m})=\frac{1}{m}\sum_{i=1}^{m}\ln\frac{Z\left(\bm{o}^{(i)}\right)}{Z}.$
(25)

To calculate $\nabla_{\bm{\Theta}}\mathcal{L}$ let us first compute
$\nabla_{\bm{\Theta}}\ln Z(\bm{o})$. Working out the derivative of the first
term

$\nabla_{\bm{\Theta}}\ln
Z(\bm{o})=\frac{1}{Z(\bm{o})}\sum_{\bm{h}}\int\mathrm{d}\bm{x}\left[-\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\right]e^{-E(\bm{x},\bm{h})}\chi(\bm{x},\bm{o}).$
(26)

Substituting Eq. (24) with
$p(\bm{o})\equiv\sum_{\bm{h}}\int\mathrm{d}\bm{x}p(\bm{o},\bm{x},\bm{h})$ in
Eq. (26)

$\frac{e^{-E(\bm{x},\bm{h})}\chi(\bm{x},\bm{o})}{Z(\bm{o})}=\frac{p(\bm{o},\bm{x},\bm{h})}{p(\bm{o})}=p(\bm{x},\bm{h}|\bm{o}),$
(27)

shows that the normalisation constant $\Xi$ of Eq. (24) cancels out exactly,
and therefore

$\nabla_{\bm{\Theta}}\ln
Z(\bm{o})=-\langle\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\rangle_{p(\bm{x},\bm{h}|\bm{o})}\,.$
(28)

For the data-independent term
$Z=\sum_{\bm{h}}\int\mathrm{d}\bm{x}\exp[-E(\bm{x},\bm{h})]$, we obtain the
standard result [24]

$\nabla_{\bm{\Theta}}\ln
Z=-\langle\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\rangle_{p(\bm{x},\bm{h})}.$
(29)

Putting the two terms together, we arrive at the desired result

$\nabla_{\bm{\Theta}}\mathcal{L}=-\left(\frac{1}{m}\sum_{i=1}^{m}\left\langle\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\right\rangle_{p(\bm{x},\bm{h}|\bm{o}^{(i)})}-\langle\nabla_{\bm{\Theta}}E(\bm{x},\bm{h})\rangle_{p(\bm{x},\bm{h})}\right).$
(30)

##  Appendix C Derivation harmonium as a mixture of Gaussians

On a high level, the event time density has four temporal modes located at
$\bm{v}^{(1)},\dots,\bm{v}^{(4)}$: two corresponding to binary colour
$x^{A}_{1}=0$ ($\bm{v}^{(1)}$ and $\bm{v}^{(2)}$) and two for colour
$x^{A}_{1}=1$ ($\bm{v}^{(3)}$ and $\bm{v}^{(4)}$). Since there are no
numerical factors $\bm{x}^{C}$ we disregard corresponding terms in $E$, so
that our goal will be to compute
$p(\bm{x}^{A},\bm{x}^{B})=\sum_{\bm{h}}\exp[-E(\bm{x}^{A},\bm{x}^{B},\bm{h})]/Z$
with weights that fit the distribution. We therefore allocate one hidden unit
$h_{i}$ for each mode $\bm{v}^{(i)}$. For convenience, write $x^{A}\equiv
x^{A}_{1}$ since there is only one binary variable (colour). Simplifying, by
setting the visible biases to zero ($\bm{c}=\bm{a}^{A}=\bm{a}^{B}=0$), we can
evaluate $p(x^{A},\bm{x}^{B}$) up to a normalisation constant $Z$ by
marginalising out $\bm{h}$:

$p(x^{A},\bm{x}^{B})=\sum_{\bm{h}}p(x^{A},\bm{x}^{B},\bm{h})\propto\sum_{\bm{h}}\exp[-E(x^{A},\bm{x}^{B},\bm{h})]=\prod_{j=1}^{4}1+e^{-\phi_{j}(x^{A},\bm{x}^{B})},$
(31)

where $\phi_{j}(x^{A},\bm{x}^{B})$ [Eq. (12)] groups energy terms proportional
to latent state $h_{j}$. Simplifying further, we substitute receptive fields
$\bm{V}$ and $\bm{W}^{B}$ in terms of its shape $\bm{\alpha}=|\bm{V}|+1$ and
rate $\bm{\beta}=\bm{W}^{B}$, and replace ${W}^{A}_{1j}$ by $w^{A}_{j}$ for
notational convenience. In this notation, we have

$\phi_{j}=-\ln(\bm{x}^{B})\cdot(\bm{\alpha}_{j}-1)+\bm{x}^{B}\cdot\bm{\beta}_{j}+w^{A}_{j}x^{A}+b_{j},$
(32)

where we used $\bm{\alpha}_{j}$ to denote column $j$ of matrix $\bm{\alpha}$,
and similarly for $\bm{\beta}$. To pin $x^{A}$ to its corresponding value
$\tilde{x}^{A}_{j}$ of mode $j$, let
$w^{A}_{j}=-q(\tilde{x}^{A}_{j}-\frac{1}{2})$ and recall the Le Roux-Bengio
Kronecker delta identity [41]:

$\lim_{q\rightarrow\infty}\exp\left[-(w^{A}_{j}x^{A}-w^{A}_{j}\tilde{x}^{A}_{j})\right]=\delta_{x^{A},\tilde{x}^{A}_{j}},$
(33)

when both $x^{A}$ and $\tilde{x}^{A}_{j}$ are binary valued. Next, observe
that most of the temporal mode’s weight are concentrated around its maximum
$v^{(j)}_{i}=(\alpha_{ij}-1)/\beta_{ij}$, justifying a Taylor expansion around
it:

$(\alpha_{ij}-1)\ln
x^{B}_{i}-\beta_{ij}x^{B}_{i}\approx(\alpha_{ij}-1)\left[\ln
v_{i}^{(j)}-1\right]-\frac{\beta_{ij}}{v^{(j)}_{i}}\frac{(x^{B}_{i}-v_{i}^{(j)})^{2}}{2}+\mathcal{O}[(x^{B}_{i}-v_{i}^{(j)})^{3}].$
(34)

With both identities [Eqs. (33) and (34)] in hand, sweep the constants in the
bias term:

$b_{j}=-w_{j}^{A}\tilde{x}^{A}_{j}+\sum_{i=1}^{2}(\alpha_{ij}-1)\left[\ln
v_{i}^{(j)}-1\right]+\Lambda,$ (35)

together with a convergence factor $\Lambda$. Substituting Eqs. (34) and (35)
in (32), and using identity (33) we have:

$e^{-\phi_{j}}\approx\delta_{x^{A},\tilde{x}^{A}_{j}}e^{-\Lambda}\exp\left[-\frac{1}{2}(\bm{x}^{B}-\bm{v}^{(j)})^{T}\bm{\Sigma}_{j}^{-1}(\bm{x}^{B}-\bm{v}^{(j)})\right],$
(36)

with
$\Sigma_{j}=\mathrm{diag}\left[\frac{\alpha_{1j}-1}{\beta_{1j}^{2}},\frac{\alpha_{2j}-1}{\beta_{2j}^{2}}\right]$
a diagonal covariance matrix. Assuming that there is little overlap between
the modes, i.e., $e^{-\phi_{j}}e^{-\phi_{k}}\approx 0$ for $j\neq k$, we have

$\tilde{p}(x^{A},\bm{x}^{B})=\prod_{j=1}^{4}1+e^{-\phi_{j}}\approx
1+\sum_{j=1}^{4}e^{-\phi_{j}}=1+2\pi
e^{-\Lambda}\sum_{j=1}^{4}\delta_{x^{A},\tilde{x}^{A}_{j}}\sqrt{|\bm{\Sigma}_{j}|}\mathcal{N}\left(\bm{v}^{(j)},\bm{\Sigma}_{j}\right),$
(37)

where $|\bm{\Sigma}_{j}|$ is used to denote the determinant of covariance
matrix $\bm{\Sigma}_{j}$. Finally, assume that each Gaussian
$\mathcal{N}\left(\bm{v}^{(j)},\bm{\Sigma}_{j}\right)$ is sufficiently
localised on the $[0,1]\times[0,1]$ unit square so that

$\int_{0}^{1}\mathrm{d}x^{B}_{1}\int_{0}^{1}\mathrm{d}x^{B}_{2}\mathcal{N}\left(\bm{v}^{(j)},\bm{\Sigma}_{j}\right)\approx\int_{-\infty}^{\infty}\mathrm{d}x^{B}_{1}\int_{-\infty}^{\infty}\mathrm{d}x^{B}_{2}\mathcal{N}\left(\bm{v}^{(j)},\bm{\Sigma}_{j}\right)=1.$
(38)

This integral identity allows us to normalise $\tilde{p}(x^{A},\bm{x}^{B})$

$\sum_{x^{A}\in\\{0,1\\}}\int_{0}^{1}\mathrm{d}x^{B}_{1}\int_{0}^{1}\mathrm{d}x^{B}_{2}\tilde{p}(x^{A},\bm{x}^{B})\approx
2+2\pi e^{-\Lambda}\sum_{j=1}^{4}\sqrt{|\bm{\Sigma}_{j}|},$ (39)

giving rise to the overall solution as a mixture of Gaussians:

$p(x^{A},\bm{x}^{B})\approx\sum_{j=1}^{4}\pi_{j}\delta_{x^{A},\tilde{x}^{A}_{j}}\mathcal{N}\left(\bm{v}^{(j)},\bm{\Sigma}_{j}\right),$
(40)

with weights
$\pi_{j}=\frac{\sqrt{|\bm{\Sigma}_{j}|}}{\sum_{k=1}^{4}\sqrt{|\bm{\Sigma}_{k}|}}$
after choosing a sufficiently large negative convergence factor $\Lambda$.
Finally, substituting $\bm{W}^{B}$ and $\bm{V}$ back into $\bm{v}^{(j)}$ and
$\bm{\Sigma}_{j}$ we arrive at the mean and the (diagonal) covariance matrix
in terms of the receptive fields:

$v_{i}^{(j)}=\frac{\|V_{ij}\|}{W^{B}_{ij}},\quad(\Sigma_{ii})_{j}=\frac{\|V_{ij}\|}{(W^{B}_{ij})^{2}}.$
(41)

##  Appendix D Experimental aspects

###  D.1 Metrics

The concordance index [1] orders the data according to the event time, and
measures the amount of data pairs in which the model’s risk prediction is
ordered concordantly. The concordance index is thus independent of the exact
risk scores but only measures their relative ranking. We therefore chose a
fixed time point $t$, and defined the risk score as the predicted survival at
that time point i.e.,

$r_{i}=S\left(x_{i}=t|\bm{o}_{-i}\right)=\frac{p\left(x_{i}>t,\bm{o}_{-i}\right)}{p(\bm{o}_{-i})},$
(42)

where $\bm{o}_{-i}$ denotes observation $\bm{o}$ with element $i$ removed. In
addition, we used $r_{i}$ to compute the Brier loss [39] to measure the
calibration at time point $t$. Notice that Eq. (42) factors in the survival
information from all other survival variables, when there is more than one
time-to-event variable. The risk score $r_{i}$ when marginalised over survival
variable $j$ is obtained by censoring at time zero [i.e.,
$o_{j}=(\xi_{j}=0,e_{j}=0)$] so that $r_{i}=S(x_{i}=t|\bm{o}_{-\\{i,j\\}})$.
Observe, moreover, that the right hand side of Eq. (42) can be evaluated in
terms of its unnormalised probabilities since the partition function cancels
out.

There are atleast two ways to compute $r_{i}$ via the unnormalised probability
density
$\tilde{p}(\bm{o},\bm{x},\bm{h})=e^{-E(\bm{x},\bm{h})}\chi(\bm{x},\bm{o})$:
(i) integrate out $\bm{x}$ analytically and then sum over $\bm{h}$ numerically
or (ii) carry out the $\bm{h}$ sum analytically and numerically marginalise
over $\bm{x}$. While the computational complexity of the former method is
linear in the number of visible units $n_{v}$ and exponential in the number of
latent states $n_{h}$, the latter scales linearly in $n_{h}$ and roughly
exponentially in the number variables with censored/missing values. We
therefore used, for the datasets presented here, method (i) when $n_{h}<10$
and method (ii) otherwise.

###  D.2 Generation of samples

To sample from Eqs. (10), (11), and (17) requires samples from the sigmoid
function, Gaussian distribution, right truncated Gamma distribution and the
interval truncated Gamma distribution. Gaussian samples can be generated using
the SciPy routine and binary states can be sampled by picking 1 when the
sigmoid activation function exceeds a $[0,1]$ uniformly sampled threshold, and
0 otherwise. To sample from the $[t_{<},1]$ interval truncated Gamma
distribution

$p_{[t_{<},1]}^{\Gamma}(x|\alpha,\beta)=\frac{x^{\alpha-1}e^{-\beta
x}}{\int^{1}_{t_{<}}\mathrm{d}t\,t^{\alpha-1}e^{-\beta
t}}=\frac{\theta(x-t_{<})}{1-t_{<}^{\alpha}\frac{\gamma^{*}(\alpha,t_{<}\beta)}{\gamma^{*}(\alpha,\beta)}}p_{\Gamma}\left(x|\alpha,\beta\right),$
(43)

and the right truncated Gamma distributions
$p_{\Gamma}\left(x|\alpha,\beta\right)$ [Eq. (21)], observe that the samples
from the latter can be obtained from the former with $t_{<}=0$. Unfortunately,
we are not aware of any existing algorithms that can generate samples from
intervals for $\alpha>0$ and both $\beta\geq 0$ and $\beta<0$. By noticing
that

$x^{\alpha-1}e^{-\beta x}\leq e^{(\alpha-1-\beta)x}\frac{1}{\exp(\alpha-1)},$
(44)

for $\alpha>1$, we propose the following rejection algorithm $\forall\beta$
and $\alpha>1$:

1:while True do

2: sample $u\sim U_{[0,1]}$ and $y\sim c_{[0,1-t_{<}]}(y|\alpha-1-\beta)$

3: compute $p_{\mathrm{accept}}(x)$ with $x=1-y$.

4: if $u\leq p_{\mathrm{accept}}$ then

5: return $x$

6: end if

7:end while

Algorithm 2 Sampling method for the $[t_{<},1]$ interval truncated gamma
distribution for $\alpha>1$.

with

$p_{\mathrm{accept}}(x)=\left(\frac{x}{\exp(x-1)}\right)^{\alpha-1},$ (45)

$U_{[0,1]}$ the uniform distribution on the unit interval, and
$c_{[0,t]}(x|\lambda)$ the exponential decaying distribution

$c_{[0,t]}(x|\lambda)=\frac{\lambda e^{-\lambda x}}{1-e^{-t\lambda}},$ (46)

normalised on the interval $[0,t]$. Samples from Eq. (46) can be generated by
inverting its cumulative distribution
$\mathcal{C}(x|\lambda,t)=(1-\exp[-\lambda x])/(1-\exp[-t\lambda])$.

###  D.3 Initialisation of parameters

For the categorical weights (in $E_{A}$), $\bm{W}^{A}$ was drawn from the
Gaussian $\mathcal{N}(0,0.01)$ and the bias $a^{A}_{i}$ was initialised to
$\ln[(1-p_{i})/p_{i}]$ where $p_{i}$ is the corresponding average value of the
categorical variable in the training set (ignoring any missing values), as
suggested in Ref. [27]. Initialisation of the parameters in $E_{C}$ followed
Ref. [29] by using Glorot-Bengio samples [42]
$\bm{W}^{C}\sim\left[-\sqrt{\frac{6}{|H|+|C|}},\sqrt{\frac{6}{|H|+|C|}}\right]$
and $\bm{\sigma}$ was treated as an adjustable parameter with initial value 1.
The bias $\bm{a}^{C}$ was set to zero since the input features can be
standardised prior to training. Similarly, for the time-to-event parameters
(in $E_{B}$) we used
$\bm{W}^{B}\sim\left[-\sqrt{\frac{6}{|H|+|B|}},\sqrt{\frac{6}{|H|+|B|}}\right]$
and sampled both $\bm{V}$ and $\bm{c}$ uniformly from
$\left[0,2\sqrt{\frac{6}{|H|+|B|}}\right]$ to ensure unit variance [42], and
picked $\bm{a}^{B}=0$. Finally, hidden biases were set to $\bm{b}=0$ as
recommended in [27].

###  D.4 Settings of benchmark real world datasets

The concordance index [1] and the Brier loss [39] (which was not available for
the SVM) were measured using 5x5 nested cross-validation [43] where the inner
loop was used to hyperparameter tune the model with the random search
algorithm [44] from Scikit-learn [45] using 50 samples. Since both nvalt
datasets consists of two time-to-event variables while the Cox model, SVM, and
random forest can only consider a single time-to-event variable, we chose to
train and evaluate these models on the OS while the harmonium was trained on
both survival variables and was evaluated on the OS. The Brier loss was
computed at $\frac{1}{2}\tau_{\mathrm{max}}$ with $\tau_{\mathrm{max}}$ the
maximum survival in the entire dataset. For the harmonium, the same time point
were used for computing risk scores (see D.1), and we set time horizon to
$\bm{\tau}_{\mathrm{max}}$. Hyperparameters were tuned to optimise the
concordance index, where the harmonium factored in the RFS and SBMFS variable
to predict OS in the nvalt8 and nvalt11 dataset, respectively.

For the Cox model from lifelines [32], the regularisation term $R(\bm{\beta})$
is parametrised as

$R=\frac{\lambda_{C}}{2}\left[(1-\ell_{1})\|\beta\|^{2}+\ell_{1}\|\beta\|\right],$
(47)

where $\bm{\beta}$ are the coefficients of the model. Parameters $\lambda_{C}$
and $\ell_{1}$ were sampled log-uniformly from the intervals
$[10^{-5},10^{3}]$ and $[10^{-5},1]$, respectively.

For the survival SVM, the hyperparameter of the squared Hinge loss $\alpha$
was sampled log uniformly from $[2^{-12},2^{12}]$ while the ranking ratio $r$
was uniformly sampled from $[0,1]$ in steps of 0.05, as suggested in Ref.
[17].

For the random survival forest, we selected a maximum tree depth of 7 (instead
of unbounded) to reduce the memory footprint, and varied (i) the number of
estimators as $2^{j}$ uniformly from $j=0,\dots,10$, (ii) the minimum of
samples required for a split as $2^{k}$ uniformly from $k=1,\dots,5$, (iii)
the minimum number of samples per leaf as $2^{l}$ uniformly from $l=0,\dots,5$
and, (iv) maximum number of features to consider per split by randomly
selecting any of $\sqrt{n}$, $\log_{2}n$, or $n$ with equal probability, with
$n$ the number of features.

Finally, for the harmonium (i) the number of hidden units, (ii) learning rate,
(iii) number of epochs to train, (iv) mini batch size, and (v) $L_{2}$ penalty
$R(\bm{\Theta})=\lambda_{H}/2\Theta^{2}$ were all sampled log uniformly from
$[1,128]$, $[10^{-5},5\cdot 10^{-2}]$, $[500,10^{5}]$, $[25,10^{3}]$, and
$[10^{-5},10^{-1}]$, respectively. In each gradient step, a part of the
previous update was retained using a momentum fraction $1-f$, where $f$ was
chosen uniformly from $[0,0.9]$. And lastly, we allowed the Gibbs chain of the
negative phase to persist [31] instead of re-initialising it each step as in
Algorithm 1. We considered this as an hyperparameter as well, and chose either
options with 50 % chance, and fixed the number of contrastive divergence steps
to 1. The following exceptions were made to these settings: (i) for the
democracy and ncctg dataset the number of epochs was capped to $5\cdot 10^{4}$
to reduce the computation time, (ii) and we lowered the maximum learning rate
for the nvalt8 dataset to 0.0125 to prevent numerical instability.

## References

  * [1] F. E. Harrell Jr, K. L. Lee, D. B. Mark, Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors, Statistics in Medicine 15 (4) (1996) 361–387.  doi:https://doi.org/10.1002/(SICI)1097-0258(19960229)15:4<361::AID-SIM168>3.0.CO;2-4.
  * [2] D. G. Kleinbaum, Survival analysis, Springer, 1996.
  * [3] G. Schwarzer, W. Vach, M. Schumacher, On the misuses of artificial neural networks for prognostic and diagnostic classification in oncology, Statistics in Medicine 19 (4) (2000) 541–561.
  * [4] T. G. Clark, M. J. Bradburn, S. B. Love, D. G. Altman, Survival analysis part i: basic concepts and first analyses, British journal of cancer 89 (2) (2003) 232–238.
  * [5] D. R. Cox, Regression models and life-tables, Journal of the Royal Statistical Society: Series B (Methodological) 34 (2) (1972) 187–202.
  * [6] H. Binder, M. Schumacher, Allowing for mandatory covariates in boosting estimation of sparse high-dimensional survival models, BMC bioinformatics 9 (1) (2008) 14.
  * [7] D. Faraggi, R. Simon, A neural network model for survival data, Statistics in Medicine 14 (1) (1995) 73–82.
  * [8] J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, Y. Kluger, Deepsurv: personalized treatment recommender system using a cox proportional hazards deep neural network, BMC medical research methodology 18 (1) (2018) 24\.
  * [9] T. Ching, X. Zhu, L. X. Garmire, Cox-nnet: an artificial neural network method for prognosis prediction of high-throughput omics data, PLoS computational biology 14 (4) (2018) e1006076.
  * [10] E. Biganzoli, P. Boracchi, L. Mariani, E. Marubini, Feed forward neural networks for the analysis of censored survival data: a partial logistic regression approach, Statistics in Medicine 17 (10) (1998) 1169–1186.
  * [11] E. Giunchiglia, A. Nemchenko, M. van der Schaar, Rnn-surv: A deep recurrent model for survival analysis, in: International Conference on Artificial Neural Networks, Springer, 2018, pp. 23–32.
  * [12] M. F. Gensheimer, B. Narasimhan, A scalable discrete-time survival model for neural networks, PeerJ 7 (2019) e6257.
  * [13] H. Ishwaran, U. B. Kogalur, E. H. Blackstone, M. S. Lauer, Random survival forests, Annals of Applied Statistics 2 (3) (2008) 841–860.
  * [14] P. K. Shivaswamy, W. Chu, M. Jansche, A support vector approach to censored targets, in: Seventh IEEE International Conference on Data Mining (ICDM 2007), IEEE, 2007, pp. 655–660.
  * [15] V. Van Belle, K. Pelckmans, J. A. K. Suykens, S. Van Huffel, Support vector machines for survival analysis, in: Proceedings of the Third International Conference on Computational Intelligence in Medicine and Healthcare (CIMED2007), 2007, pp. 1–8.
  * [16] V. Van Belle, K. Pelckmans, S. Van Huffel, J. A. K. Suykens, Support vector methods for survival analysis: a comparison between ranking and regression approaches, Artificial intelligence in medicine 53 (2) (2011) 107–118.
  * [17] S. Pölsterl, N. Navab, A. Katouzian, Fast training of support vector machines for survival analysis, in: Joint European Conference on Machine Learning and Knowledge Discovery in Databases, Springer, 2015, pp. 243–259.
  * [18] L.-J. Wei, D. Y. Lin, L. Weissfeld, Regression analysis of multivariate incomplete failure time data by modeling marginal distributions, Journal of the American statistical association 84 (408) (1989) 1065–1073.
  * [19] L. J. Wei, , D. V. Glidden, An overview of statistical methods for multiple failure time data in clinical trials, Statistics in Medicine 16 (8) (1997) 833–839.
  * [20] D. O. Dean, D. J. Bauer, M. J. Shanahan, A discrete-time multiple event process survival mixture (mepsum) model., Psychological Methods 19 (2) (2014) 251.
  * [21] W. Schnedler, Likelihood estimation for censored random vectors, Econometric Reviews 24 (2) (2005) 195–217.
  * [22] P. Smolensky, Information processing in dynamical systems: Foundations of harmony theory, Tech. rep., Colorado Univ at Boulder Dept of Computer Science (1986).
  * [23] G. E. Hinton, Training products of experts by minimizing contrastive divergence, Neural computation 14 (8) (2002) 1771–1800.
  * [24] I. Goodfellow, Y. Bengio, A. Courville, Deep learning, MIT press, 2016.
  * [25] G. E. Hinton, T. J. Sejnowski, Learning and relearning in boltzmann machines, Parallel distributed processing: Explorations in the microstructure of cognition 1 (282-317) (1986) 2.
  * [26] P. M. Long, R. A. Servedio, Restricted boltzmann machines are hard to approximately evaluate or simulate, in: Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, Omnipress, Madison, WI, USA, 2010, p. 703–710.
  * [27] G. E. Hinton, A practical guide to training restricted boltzmann machines, in: Neural networks: Tricks of the trade, Springer, 2012, pp. 599–619.
  * [28] L. Theis, S. Gerwinn, F. Sinz, M. Bethge, In all likelihood, deep belief is not enough, The Journal of Machine Learning Research 12 (Nov) (2011) 3071–3096.
  * [29] J. Melchior, N. Wang, L. Wiskott, Gaussian-binary restricted boltzmann machines for modeling natural image statistics, PloS one 12 (2) (2017).
  * [30] T. Tran, D. Phung, S. Venkatesh, Thurstonian boltzmann machines: learning from multiple inequalities, in: International Conference on Machine Learning, 2013, pp. 46–54.
  * [31] T. Tieleman, Training restricted boltzmann machines using approximations to the likelihood gradient, in: Proceedings of the 25th international conference on Machine learning, ACM, 2008, pp. 1064–1071.  doi:10.1145/1390156.1390290.
  * [32] C. Davidson-Pilon, J. Kalderstam, N. Jacobson, sean reed, B. Kuhn, P. Zivich, M. Williamson, AbdealiJK, D. Datta, A. Fiore-Gartland, A. Parij, D. WIlson, Gabriel, L. Moneda, A. Moncada-Torres, K. Stark, H. Gadgil, Jona, K. Singaravelan, L. Besson, M. S. Peña, S. Anton, A. Klintberg, GrowthJeff, J. Noorbakhsh, M. Begun, R. Kumar, S. Hussey, D. Golland, jlim13, lifelines: v0.25.5 (Sep. 2020).  doi:10.5281/zenodo.4050560.
URL https://doi.org/10.5281/zenodo.4050560

  * [33] P. H. Rossi, R. A. Berk, K. J. Lenihan, Money, Work and Crime: Some Experimental Results, New York: Academic Press, 1980.
  * [34] J. A. Cheibub, J. Gandhi, J. R. Vreeland, Democracy and dictatorship revisited, Public choice 143 (1-2) (2010) 67–101.
  * [35] M. Schumacher, G. Bastert, H. Bojar, K. Huebner, M. Olschewski, W. Sauerbrei, C. Schmoor, C. Beyerle, R. L. Neumann, H. F. Rauschecker, Randomized 2 x 2 trial evaluating hormonal treatment and the duration of chemotherapy in node-positive breast cancer patients. german breast cancer study group., Journal of Clinical Oncology 12 (10) (1994) 2086–2093.
  * [36] C. L. Loprinzi, J. A. Laurie, H. S. Wieand, J. E. Krook, P. J. Novotny, J. W. Kugler, J. Bartel, M. Law, M. Bateman, N. E. Klatt, Prospective evaluation of prognostic variables from patient-completed questionnaires. north central cancer treatment group., Journal of Clinical Oncology 12 (3) (1994) 601–607.
  * [37] D. De Ruysscher, A. M. C. Dingemans, J. Praag, J. Belderbos, C. Tissing-Tan, J. Herder, T. Haitjema, F. Ubbels, F. Lagerwaard, S. Y. El Sharouni, J. A. Stigt, E. Smit, H. van Tinteren, V. van der Noort, H. J. M. Groen, Prophylactic cranial irradiation versus observation in radically treated stage iii non–small-cell lung cancer: A randomized phase iii nvalt-11/dlcrg-02 study, Journal of Clinical Oncology 36 (23) (2018) 2366–2377.
  * [38] H. J. M. Groen, E. H. F. M. van der Heijden, T. J. Klinkenberg, B. Biesma, J. Aerts, A. Verhagen, C. Kloosterziel, R. Pieterman, B. van den Borne, H. J. M. Smit, O. Hoekstra, F. M. N. H. Schramel, V. van der Noort, H. van Tinteren, E. F. Smit, A.-M. C. Dingemans, Randomised phase 3 study of adjuvant chemotherapy with or without nadroparin in patients with completely resected non-small-cell lung cancer: the nvalt-8 study, British journal of cancer 121 (5) (2019) 372–377.
  * [39] E. Graf, C. Schmoor, W. Sauerbrei, M. Schumacher, Assessment and comparison of prognostic classification schemes for survival data, Statistics in Medicine 18 (17-18) (1999) 2529–2545.
  * [40] R. B. Paris, Incomplete gamma and related functions, NIST digital library of mathematical functions (2010) 173–192.
  * [41] N. Le Roux, Y. Bengio, Representational power of restricted boltzmann machines and deep belief networks, Neural computation 20 (6) (2008) 1631–1649.
  * [42] X. Glorot, Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in: Proceedings of the thirteenth international conference on artificial intelligence and statistics, 2010, pp. 249–256.
  * [43] M. Stone, Cross-validatory choice and assessment of statistical predictions, Journal of the Royal Statistical Society: Series B (Methodological) 36 (2) (1974) 111–133.
  * [44] J. Bergstra, Y. Bengio, Random search for hyper-parameter optimization, Journal of Machine Learning Research 13 (1) (2012) 281–305.
  * [45] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: Machine learning in Python, The Journal of Machine Learning Research 12 (2011) 2825–2830.