# Latent variable model for high-dimensional point process
with structured missingness

Maksim Sinelnikov  Manuel Haussmann  Harri Lähdesmäki

###### Abstract

Longitudinal data are important in numerous fields, such as healthcare,
sociology and seismology, but real-world datasets present notable challenges
for practitioners because they can be high-dimensional, contain structured
missingness patterns, and measurement time points can be governed by an
unknown stochastic process. While various solutions have been suggested, the
majority of them have been designed to account for only one of these
challenges. In this work, we propose a flexible and efficient latent-variable
model that is capable of addressing all these limitations. Our approach
utilizes Gaussian processes to capture temporal correlations between samples
and their associated missingness masks as well as to model the underlying
point process. We construct our model as a variational autoencoder together
with deep neural network parameterised encoder and decoder models, and develop
a scalable amortised variational inference approach for efficient model
training. We demonstrate competitive performance using both simulated and real
datasets.

##  1 Introduction

Longitudinal data arise in many domains such as healthcare, sociology and
seismology (Liu, 2015). These datasets consist of repeated measurements of
unique instances, e.g. patients, collected over time. However, real-world
applications pose several challenges for practitioners: measurements are
typically high-dimensional and contain non-trivial missingness patterns, and
time points of the observations are not deterministic but rather arise from an
unknown stochastic process. These challenges are characteristic of many real
biomedical datasets, such as electronic health records.

Variational autoencoders (VAEs) have become a popular approach to model high-
dimensional data (Kingma & Welling, 2014; Rezende et al., 2014). A notable
limitation of standard VAEs is their assumption that latent variables
factorize across samples, hence ignoring correlations between observations and
making the models inappropriate for temporal and longitudinal datasets.
Several recent works (Casale et al., 2018; Fortuin et al., 2020; Ramchandran
et al., 2021) have addressed this issue by incorporating Gaussian process (GP)
priors for the latent variables, creating a probabilistic model that is
capable of modelling arbitrary correlations between latent encodings.

The simplest form of missingness is missing completely at random (MCAR), i.e.,
the missingness pattern is independent of the observed and unobserved data.
While it is generally feasible to handle MCAR in most latent-variable models,
more complex patterns of structured missingness (Rubin, 1976; Mitra et al.,
2023) require additional modeling capacities. Extending VAEs to be able to
model various structured missingness patterns has recently become the focus of
several papers (Collier et al., 2020; Ipsen et al., 2021; Ghalebikesabi et
al., 2021). However, these methods still lack the ability to model
correlations among observations or missingness masks, thus limiting their
applicability to temporal data.

While VAEs have been applied to various biomedical datasets, the existing
methods cannot consider observation time points as random variables. Instead,
they have to treat time as a deterministic covariate and, therefore, loose
useful information embedded in its stochastic nature. A separate line of
research has proposed methods to model unknown temporal point processes
primarily using GP-based methods (Lloyd et al., 2015; Liu & Hauskrecht, 2019).
Overall, the field lacks versatile modeling methods that would allow modeling
high-dimensional, marked point-processes that may be corrupted by structured
missingness.

#### Contributions.

In this work, we propose a novel deep latent variable model (DLVM), that is
specifically designed to capture structured missingness and uses temporal
point processes to model time. We construct the model by introducing three
sets of latent variables with GP priors, which model observations, missingness
masks and point process. To adapt the model for longitudinal data, we rely on
longitudinal additive kernels (Ramchandran et al., 2021) for the latent
representations of observations and masks. Additionally, to use the
information embedded in the temporal point process, we provide the intensity
of the point process as an additional input to the GP kernel functions of
observations and missingness masks. See Figure 1 for a high-level summary of
our model and Appendix G for an extended visualization. We present two
variations of our proposal, the simplified _longitudinal latent variable model
with structured missingness_ (LLSM) without a temporal point process, and the
full model, _longitudinal latent variable model for high-dimensional point
process with structured missingness_ (LLPPSM). To summarize, our contributions
are that

  * _(i)_

we present a latent variable model able to capture structured missingness in
the context of longitudinal data;

  * _(ii)_

we extend this model by a temporal point process and use the inferred
intensity of the process as an additional input to the model;

  * _(iii)_

we compare the performance of our two model variants against baseline methods
on several datasets and report state-of-the-art results on a variety of tasks.

Figure 1: Conceptual overview of our model. Shaded, partially-shaded, and
blank rectangles refer to observed, partially observed, and latent components.

##  2 Related Work

We summarize and compare previous methods in Table 1.

Table 1: A summary of related methods.

Model Temporal structure Other covariates Structured missingness Modelling
timestamps Generative Reference VAE missing ✗ ✗ ✓ ✗ ✓ Collier et al. (2020)
not-MIWAE ✗ ✗ ✓ ✗ ✓ Ipsen et al. (2021) GRUI-GAN ✓ ✗ ✓ ✗ ✓ Luo et al. (2018)
GPPVAE ✓ Limited ✗ ✗ ✓ Casale et al. (2018) GP-VAE ✓ ✗ ✗ ✗ ✓ Fortuin et al.
(2020) L-VAE ✓ ✓ ✗ ✗ ✓ Ramchandran et al. (2021) GPRPP ✓ ✗ ✗ ✓ Limited to
timestamps Liu & Hauskrecht (2019) LLSM ✓ ✓ ✓ ✗ ✓ This work LLPPSM ✓ ✓ ✓ ✓ ✓
This work

#### Challenges with missing data.

In the pioneer work, Rubin (1976) identified three classes of missing data:
missing completely at random (MCAR), missing at random (MAR), and missing not
at random (MNAR). For MCAR, the missingness mechanism is independent of both
observed and unobserved variables. In case of MAR, the missingness depends on
the observed attributes. Whereas if data is MNAR, missing readings are
dependent on the unobserved data, or systematic factors that are not accounted
for in the experiment. The last two are examples of structured missingness.
Although MCAR case can be handled by simply excluding missing elements from
the analysis without introducing a bias, the same does not hold for the other
scenarios. Despite the undoubted utility of Rubin’s taxonomies, Mitra et al.
(2023) emphasized that they do not fully account for high-dimensional patterns
of structured missingness, frequently encountered in modern ML applications.
They also proposed a set of current grand challenges in learning from data
with structured missingness and claimed that the field of learning with
missing values needs to be further advanced.

#### DLVMs for missing data.

Various methods have been proposed to address the challenge of missing values
within generative models. Collier et al. (2020) employed variational
autoencoder by concatenating input with the missing mask. While this approach
allows to model MAR and MNAR scenarios, it fails to model temporal
correlations and does not take into account auxiliary covariate information.
Mattei & Frellsen (2019) used an idea of importance sampling to derive a
Missing Importance Weighted Autoencoder(MIWAE) bound for training DLVMs under
MAR condition. Building upon this work, Ipsen et al. (2021) expanded this
method to MNAR scenario by directly modeling the missing mask. However, both
these approaches lack the ability to model temporal correlation in the latent
space, hence are not suitable for longitudinal setting.

#### Temporal DLVMs.

To model temporal data, a set of methodologies has emerged, exploring the use
of GP priors for latent variables. Casale et al. (2018) introduced the GPPVAE
model to integrate both view and object information in a GP prior through the
product kernel. Fortuin et al. (2020) proposed the GP-VAE that assigns
individual GP prior for the time-series of each subject. While these methods
allow to model subject-specific temporal structure, they have limited or no
functionality to account for possible other auxiliary covariates. Ramchandran
et al. (2021) introduced L-VAE, a model that uses a multi-output additive GP
prior and is well-suited for longitudinal data by leveraging carefully
designed interaction kernels. The main drawback of all these approaches is
their limitation to MCAR modeling which is naïve in many domains, such as
healthcare.

Another direction of research to deal with temporal data focuses on recurrent
neural networks (RNNs). Che et al. (2018) developed GRU-D that incorporates
_masking_ and _time interval_ into a deep model architecture, making it
possible to model structured missingness patterns. They also proposed to use
decaying mechanism to handle irregularly-sampled timestamps. Luo et al. (2018)
extended this model for time-series imputation by employing generative
adversarial networks (GRUI-GAN). However, it is not straightforward to
incorporate auxiliary information into these models.

#### Temporal Point Process.

Modelling temporal point processes has been a subject of several studies in
recent years. Classical statistical approaches (Puri & Tuan, 1986) use maximum
likelihood estimation to infer parameters of a model. However, they require
specifying a parametric form of the intensity function which significantly
limits their applications. Neural network based models typically employ RNNs
(Du et al., 2016). Lloyd et al. (2015) proposed to model intensity function
with Gaussian processes and squared link function that leads to closed-form
solution. John & Hensman (2018) extended this model to variational Fourier
features. Both approaches use the current time point, not taking into account
the previous history of events which limits them to model only inhomogeneous
Poisson processes, a potentially unrealistic assumption in real-world
scenarios. Liu & Hauskrecht (2019) overcome this problem by incorporating the
previous $D$ timestamps into the computation of GP kernel function.

##  3 Methods

###  3.1 Background

#### Problem setup.

We are given $N=\sum_{p=1}^{P}n_{p}$ observations, where $P$ denotes the
number of unique instances, e.g., patients, and $n_{p}$ is the number of
observations of instance $p$. The longitudinal _response variables_ (or
_marks_) of instance $p$ are denoted as
$\boldsymbol{y}_{p}=[y_{1}^{p},\ldots,y_{n_{p}}^{p}]\in\mathbb{R}^{K\times
n_{p}}$, where each sample $y_{i}^{p}\in\mathcal{Y}=\mathbb{R}^{K}$. Each
sample $y_{i}\in\mathcal{Y}$ can be split into observed and missing parts,
$y_{i}^{\text{o}}$ and $y_{i}^{\text{m}}$, with a corresponding _binary mask_
$m_{i}\in\\{0,1\\}^{K}$ specifying which features of $y_{i}$ are missing (1 is
observed, 0 is missing). The _auxiliary covariate_ information of instance $p$
is $\boldsymbol{x}_{p}=[x_{1}^{p},\ldots,x_{n_{p}}^{p}]$, where each
$x_{t}^{p}\in\mathcal{X}=\mathcal{X}_{1}\times\cdots\times\mathcal{X}_{Q}$ is
$Q$-dimensional vector. Covariates can be both discrete and continuous and
represent, for instance, patient’s age, gender, etc. For notational
convenience, we separately denote the _time points_ of the measurements of
subject $p$ as $\boldsymbol{t}^{p}=[t_{1}^{p},\ldots,t_{n_{p}}^{p}]^{T}$ that
are considered as random variables, and $\mathcal{X}^{\text{static}}$ as the
set of covariates that do not depend on time, such that $t\subset\mathcal{X}$
and $\mathcal{X}^{\text{static}}\subset\mathcal{X}$. Observations from all $P$
instances (e.g. patients) form our longitudinal data matrix $\boldsymbol{y}$,
matrix of missing values $\boldsymbol{m}$, covariate matrix $\boldsymbol{x}$,
and vector of timestamps $\boldsymbol{t}$, defined respectively as

$\displaystyle\boldsymbol{y}$
$\displaystyle=[\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{P}]=[y_{1},\ldots,y_{N}]$
$\displaystyle\in\mathbb{R}^{K\times N},$ $\displaystyle\boldsymbol{m}$
$\displaystyle=[\boldsymbol{m}_{1},\ldots,\boldsymbol{m}_{P}]=[m_{1},\ldots,m_{N}]$
$\displaystyle\in\\{0,1\\}^{K\times N},$ $\displaystyle\boldsymbol{x}$
$\displaystyle=[\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{P}]=[x_{1},\ldots,x_{N}]$
$\displaystyle(\text{size }Q\times N),$ $\displaystyle\boldsymbol{t}$
$\displaystyle=[\boldsymbol{t}_{1}^{T},\ldots,\boldsymbol{t}_{P}^{T}]^{T}=[t_{1},\ldots,t_{N}]^{T}$
$\displaystyle\in\mathbb{R}^{N}.$

We rely on a latent space ${\mathcal{Z}=\mathbb{R}^{L}}$ and combine the
latent embedding for all $N$ samples as
${\boldsymbol{z}=[z_{1},\ldots,z_{N}]\in\mathbb{R}^{L\times N}}$.

#### Variational autoencoders.

Assuming a deep latent variable model
$p_{\omega}(y,z)=p_{\psi}(y|z)p_{\theta}(z)$, inference of the posterior
$p_{\omega}(z|y)=p_{\psi}(y|z)p_{\theta}(z)/p_{\omega}(y)$ is in general
intractable, as the evidence $p_{\omega}(y)$ cannot be computed analytically
due to the highly non-linear relationship between $z$ and $y$. Common practice
is to rely on amortized inference (Kingma & Welling, 2014; Rezende et al.,
2014), where a parameterized approximation, $q_{\phi}(z|y)$ to the true
posterior is inferred by optimizing a lower bound to the evidence,

$\log p_{\omega}(y)\geq{\operatorname*{\mathbb{E}}}_{q}\left[\log
p_{\psi}(y|z)\right]-\text{KL}(q_{\phi}(z|y)||p_{\theta}(z)),$

with respect to all parameters. Here, KL is the Kullback-Leibler divergence.
Usually, likelihood $p_{\psi}(y|z)$, prior $p_{\theta}(z)$, and variational
posterior $q_{\phi}(z|y)$ are assumed to be mean-field, i.e., to factorize
over their respective random variables.

#### GP-prior variational autoencoder.

Despite the computational efficiency provided by a factorized prior
$p_{\theta}(\boldsymbol{z})$ over the latent variables, its major limitation
is the inability to model correlations between data samples. Prior work
addressed this by combining VAEs with GPs, creating a powerful probabilistic
model for this task (Casale et al., 2018; Fortuin et al., 2020; Ramchandran et
al., 2021). The key difference is that the factorized prior
$p_{\theta}(\boldsymbol{z})$ is replaced by a GP-prior
$p_{\theta}(\boldsymbol{z}|\boldsymbol{x})$ which depends on auxiliary
information $\boldsymbol{x}$. The conditional generative model is then given
as

$p_{\omega}(\boldsymbol{y}|\boldsymbol{x})=\int{\prod_{i=1}^{N}{p_{\psi}(y_{i}|z_{i})}p_{\theta}(\boldsymbol{z}|\boldsymbol{x})d\boldsymbol{z}}.$

Defining a mapping from the covariates to the latent space,
$\textbf{f}:\mathcal{X}\rightarrow\mathcal{Z}$, such that
$z=\textbf{f}(x)=[f_{1}(x),\ldots,f_{L}(x)]^{T}$, these models assume a GP-
prior over each $f_{l}$,

$f_{l}(x)\sim\mathcal{GP}\left(\big{.}\mu_{l}(x),k_{l}(x,x^{\prime}|\theta_{l})\right),$

where $\mu_{l}(x)$ is the mean function and $k_{l}(x,x^{\prime}|\theta_{l})$
is the covariance function parameterized by $\theta_{l}$.

Given this prior, the $l$-th latent variable
${\bar{\boldsymbol{z}}_{l}=f_{l}(\boldsymbol{x})=[f_{l}(x_{1}),\ldots,f_{l}(x_{N})]^{T}}$
has a multivariate Gaussian distribution across the $N$ data samples

$p_{\theta_{l}}(\bar{\boldsymbol{z}}_{l}|\boldsymbol{x})=p_{\theta_{l}}(f_{l}(\boldsymbol{x}))=\mathcal{N}\Big{(}\bar{\boldsymbol{z}}_{l}\big{|}\boldsymbol{0},K^{(l)}_{\boldsymbol{x}\boldsymbol{x}}\Big{)},$

where $K^{(l)}_{\boldsymbol{x}\boldsymbol{x}}$ is a $N\times N$ covariance
matrix such that
$\\{K^{(l)}_{\boldsymbol{x}\boldsymbol{x}}\\}_{ij}=k_{l}(x_{i},x_{j}|\theta_{l})$.
We follow common practice and factorize the GP-prior across its $L$
dimensions, such that the conditional prior is given as

$p_{\theta}(\boldsymbol{z}|\boldsymbol{x})=\prod_{l=1}^{L}{p_{\theta_{l}}(\bar{\boldsymbol{z}}_{l}|\boldsymbol{x})}=\prod_{l=1}^{L}{\mathcal{N}\left(\bar{\boldsymbol{z}}_{l}\big{|}\boldsymbol{0},K^{(l)}_{\boldsymbol{x}\boldsymbol{x}}\right)}.$

The main distinction among previous GP-prior models lies in the choice of
covariance functions. GPPVAE (Casale et al., 2018) and GP-VAE (Fortuin et al.,
2020) both rely on restricted kernels that do not adequately model
longitudinal data structure. Instead, we adopt the proposal by Ramchandran et
al. (2021) who introduce a flexible additive kernel structure that is
specifically designed for longitudinal data and is capable of employing
various interactions between continuous and categorical covariates

$k_{l}(x,x^{\prime}|\theta_{l})=\sum_{r=1}^{R}{k_{l,r}(x,x^{\prime}|\theta_{l,r})}+\sigma_{zl}^{2},$

such that

$p_{\theta}(\boldsymbol{z}|\boldsymbol{x})=\prod_{l=1}^{L}{\mathcal{N}\left(\bar{\boldsymbol{z}}_{l}\Big{|}\boldsymbol{0},\sum_{r=1}^{R}{K^{(l,r)}_{\boldsymbol{x}\boldsymbol{x}}+\sigma_{zl}^{2}I_{N}}\right)}.$

###  3.2 Modeling structured missingness

These GP-prior models are capable of dealing with missing values solely by
substituting zeros or alternative values and propagating $y$ through encoder-
decoder structure to perform imputation. However, this approach lacks the
ability to model specific missingness patterns therefore making them suitable
only for an MCAR scenario, an unrealistic assumption in many real
applications. In this work, we solve this constraint and propose the
_longitudinal latent variable model with structured missingness (LLSM)_.

To model non-random missingness in VAE models various approaches exist. For
example, Mattei & Frellsen (2019) and Ipsen et al. (2021) model the dependency
between $y$ and $m$ directly, whereas Collier et al. (2020) propose a VAE
model that incorporates an additional latent variable to account for
structured missingness. In this work, we follow Collier et al. (2020) and
introduce a second latent variable $z^{m}\in\mathbb{R}^{L_{m}}$ associated
with a missingness mask $m$, while referring to the latent variables
associated to $y$ from now on as $z^{y}\in\mathbb{R}^{L_{y}}$. To properly
model MNAR we assume that $m$ depends on both $z^{m}$ and $z^{y}$. The joint
likelihood for one sample is given as

$\displaystyle p_{\omega}(y^{\text{o}},z^{y},m,z^{m}|x)$
$\displaystyle=p_{\psi_{y}}(y^{\text{o}}|z^{y},z^{m},m)p_{\psi_{m}}(m|z^{y},z^{m})$
$\displaystyle\quad\cdot p_{\theta_{y}}(z^{y}|x)p_{\theta_{m}}(z^{m}|x),$ (1)

where $y^{\text{o}}$ refers to the observed features specified by $m$. Also,
by conditioning $y$ on $z^{m}$, we can use any information, contained in the
missing mask, for instance, during an imputation task. Depending on the
application this dependency could be removed. To model correlation within the
missingness patterns, e.g., across time, or within a patient, we assign
$z^{m}$ a GP prior as well,

$z^{m}(x)\sim\mathcal{GP}\left(\big{.}0,k(x,x^{\prime}|\theta_{m})\right).$

Additionally, we assume that $y^{o}$ and $m$ are distributed as

$\displaystyle p_{\psi_{y}}(y^{\text{o}}|z^{y},z^{m},m)$
$\displaystyle=\mathcal{N}(y|g_{\psi_{y}}(z^{y},z^{m}),\boldsymbol{\sigma}^{2})\odot
m$ $\displaystyle p_{\psi_{m}}(m|z^{y},z^{m})$
$\displaystyle=\mathcal{B}\textit{er}(m|g_{\psi_{m}}(z^{y},z^{m})),$

where the decoder functions $g_{\psi_{y}}$ and $g_{\psi_{m}}$ are
parameterized by neural networks, $\odot$ denotes an element-wise Hadamard
product, and the observational variance parameters
${\boldsymbol{\sigma}^{2}=\mathrm{diag}(\sigma^{2}_{y_{1}},\ldots,\sigma^{2}_{y_{K}})}$
are optimized jointly with all other parameters via gradient descent. The
graphical model of LLSM is shown in Figure 2.

$x$$z^{m}$$z^{y}$$m$$y$ Figure 2: _The graphical model of LLSM._ Shaded,
partially-shaded, and blank circles refer to observed, partially observed, and
latent variables. The dashed arrow highlights an optional dependency.

We approximate the intractable posterior across all $N$ samples
$p(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x},\boldsymbol{y}^{\text{o}},\boldsymbol{m})$
using a mean-field amortized inference distribution

$\displaystyle
q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})$
$\displaystyle=q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})$
$\displaystyle=\prod_{i=1}^{N}{q_{\phi_{y}}(z_{i}^{y}|y^{\text{o}}_{i})q_{\phi_{m}}(z_{i}^{m}|m_{i})}$
(2)

where $q_{\phi_{y}}(z_{i}^{y}|y^{\text{o}}_{i})$ and
$q_{\phi_{m}}(z_{i}^{m}|m_{i})$ are diagonal Gaussian distributions
parameterized by neural network-based mappings from the corresponding inputs
$y^{\text{o}}_{i}$ and $m_{i}$. The lower bound on the evidence to be
optimized is given as

$\displaystyle\log p(\boldsymbol{y}^{o},\boldsymbol{m}|\boldsymbol{x})$
$\displaystyle\geq{\operatorname*{\mathbb{E}}}_{q_{y}}\left[\log{p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})}\right]$
$\displaystyle\quad+{\operatorname*{\mathbb{E}}}_{q_{m}}\left[\log{p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}\right]$
$\displaystyle\quad-\hat{\text{KL}}(q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})||p_{\theta_{y}}(\boldsymbol{z}^{y}|\boldsymbol{x}))$
$\displaystyle\quad-\hat{\text{KL}}(q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})||p_{\theta_{m}}(\boldsymbol{z}^{m}|\boldsymbol{x})),$

where, in order to maintain computational tractability, and to be able to
perform mini-batching, we substitute the KL terms with the corresponding upper
bounds $\hat{\text{KL}}$ derived by Ramchandran et al. (2021) that, for
longitudinal data, are tighter than the well-known bound proposed by Titsias
(2009). Further details on the lower bound and KL upper bound are given in
Appendix A and Appendix B.

###  3.3 Modeling Time

Prior work relying on GP-based priors suffers from a second restriction. They
often rely on time $t$ as a primary, or even the only (Fortuin et al., 2020),
covariate that is used in the covariance function. This implies that the
similarity between two measurements $y$ and $y^{\prime}$ is directly
contingent on the corresponding $t$ and $t^{\prime}$, e.g., their temporal
difference when employing a stationary kernel. While this assumption may be
reasonable in some cases, it is not universally true. For instance, consider a
scenario where a patient develops a disease, and since the progression varies
for each individual, it may be more appropriate for similarity to be
determined not solely by the elapsed time since the onset of the disease, but
rather by their current well-being, which can be captured by factors such as
the frequency of doctor visits. To properly account for such variations, we
model $t$ by a temporal point process and add the intensity of this point
process as an additional input to the GP kernel computation.

#### Temporal point processes (TPP).

A TPP (Cox & Isham, 1980) is a stochastic process over variable-length
sequences in some time interval $\mathcal{T}=\left[0,T\right]$ defined via an
intensity function $\lambda(t)$. The probability density function of $N$
observed points ${\mathcal{\boldsymbol{t}}=\\{t_{i}\in\mathcal{T}\\}}$ is then
defined as

$p(\mathcal{\boldsymbol{t}}|\lambda)=\exp\Big{(}-\int_{\mathcal{T}}{\lambda(t)dt}\Big{)}\prod_{i=1}^{N}{\lambda(t_{i})}.$

TPPs can be divided into roughly two classes: homogeneous Poisson processes,
where the intensity function only depends on the current time point $t$, and
self-exciting processes, where the occurrence of events changes the intensity.
One type of these self-exciting processes, known as the _Hawkes process_
(Hawkes, 1971), has an intensity function

$\lambda(t)=\mu+\sum_{t_{j}<t}{\nu(t-t_{j})},$

where $\nu$ is a _triggering kernel_ that characterizes the influence of past
events on intensity at time $t$ and $\mu$ is a corresponding baseline.
Inspired by the broad applicability of Hawkes processes (Hawkes, 2018), we
adopt the proposal by Liu & Hauskrecht (2019) to model such self-exciting
processes with GPs by computing kernels from the last $D$ timestamps.

#### GP point processes.

Given a latent variable $z^{\lambda}$ with

$\displaystyle\begin{split}z^{\lambda}(t)&\sim\mathcal{GP}\big{(}0,k_{\theta_{\lambda}}(v_{D},{v_{D}}^{\prime})\big{)}\\\
v_{D}&=t-t_{D},\end{split}$ (3)

where $t_{D}$ denotes $D$ previous timestamps before $t$, $v_{D}$ are the
elapsed times between $t$ and $t_{D}$, and $k_{\theta_{\lambda}}$ is an
additive kernel structure, we model the intensity as

$\lambda(t)=(z^{\lambda}(t)+\beta)^{2}$ (4)

where $\beta$ is either a trainable baseline or a function that can depend on
static covariates (John & Hensman, 2018). We choose a squared link function
due to the analytical tractability of integral that it provides (Lloyd et al.,
2015).

The posterior distribution of the intensity,

$p(\lambda|\boldsymbol{t})=\frac{p(\lambda)\exp\left(-\int_{\mathcal{T}}{\lambda(t)dt}\right)\prod_{i=1}^{N}{\lambda(t_{i})}}{\int{p(\lambda)\exp\left(-\int_{\mathcal{T}}{\lambda(t)}dt\right)\prod_{i=1}^{N}{\lambda(t_{i})}d\lambda}},$

is intractable due to the integration over $\lambda$. To overcome this
challenge, we approximate it with a variational distribution
$p(z^{\lambda}|\boldsymbol{u})q(\boldsymbol{u})$ that relies on inducing
points $\boldsymbol{u}$ for additional scalability (Quiñonero Candela &
Rasmussen, 2005). See Appendix C for a more detailed discussion.

#### LLPPSM.

Combining such a point process with our LLSM model allows us to capture
intricate missingness patterns and to effectively leverage information
embedded in the time points. We call this model _longitudinal latent variable
model for high-dimensional point process with structured missingness
(LLPPSM)_. See Figure 3 for its plate diagram.

$t_{D}$$z^{\lambda}$$\lambda$$t$$x$$x_{p}^{\text{s}}$$z^{m}$$z^{y}$$m$$y$
Figure 3: _The graphical model of LLPPSM_. (Partially) shaded, and blank
circles refer to (partially) observed, and latent variables. The dashed arrow
highlights an optional dependency, $x_{p}^{\text{s}}$ are static covariates
and $t_{D}$ the $D$ previous time steps to $t$.

Defining $z^{\lambda}$ and $\lambda$ as in Equations 3 and 4, we extend the GP
priors for $z^{y}$ and $z^{m}$ by letting their covariance kernel depend on
the rate of the TPP $\lambda$ as well, i.e.,

$z^{y}(x,\lambda)\sim\mathcal{GP}\big{(}0,k((x,\lambda(t)),(x^{\prime},\lambda(t^{\prime}))|\theta_{y})\big{)},$

and analogously for $z^{m}$. As inference of the full model remains
intractable, we once again rely on variational inference and use the following
variational approximation

$\displaystyle
q(\boldsymbol{z}^{y},\boldsymbol{z}^{m},z^{\lambda},\boldsymbol{u}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})$
$\displaystyle\qquad=q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})p(z^{\lambda}|\boldsymbol{u})q(\boldsymbol{u}),$

where $q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})$ and
$q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})$ are defined as in Equation
2, and $\boldsymbol{u}$ are the inducing points of $z^{\lambda}$. The bound to
be optimized is then given as

$\displaystyle\log
p(\boldsymbol{y}^{o},\boldsymbol{m},\boldsymbol{t}|\boldsymbol{x}^{\text{s}})\geq$
$\displaystyle\quad{\operatorname*{\mathbb{E}}}_{q}\left[\log\left(p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}\big{|}\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})p_{\psi_{m}}(\boldsymbol{m}\big{|}\boldsymbol{z}^{y},\boldsymbol{z}^{m})p(\boldsymbol{t}\big{|}\lambda)\right)\right]$
$\displaystyle\qquad-\text{KL}(q(\boldsymbol{u})||p(\boldsymbol{u}))$
$\displaystyle\qquad-{\operatorname*{\mathbb{E}}}_{p(z^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})}[\hat{\text{KL}}(q_{\phi_{y}}(\boldsymbol{z}^{y}\big{|}\boldsymbol{y}^{\text{o}})||p_{\theta_{y}}(\boldsymbol{z}^{y}\big{|}\boldsymbol{x},\lambda(\boldsymbol{t}))]$
$\displaystyle\qquad-{\operatorname*{\mathbb{E}}}_{p(z^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})}[\hat{\text{KL}}(q_{\phi_{m}}(\boldsymbol{z}^{m}\big{|}\boldsymbol{m})||p_{\theta_{m}}(\boldsymbol{z}^{m}\big{|}\boldsymbol{x},\lambda(\boldsymbol{t}))],$

where $\boldsymbol{x}^{\text{s}}$ refers to static covariates and
$\boldsymbol{x}$ is composed of $\boldsymbol{x}^{\text{s}}$ and
$\boldsymbol{t}$. See Appendix C for a detailed derivation and discussion.

###  3.4 Imputation and future prediction

Our method can be employed for various tasks such as imputation and future
prediction. The imputation is done by substituting missing elements with some
intermediate values (in our case zeros), and propagating them through the
encoder-decoder structure of our model so that the decoder imputes the missing
values.

Future predictions are obtained by evaluating the posterior predictive
distribution $p(y_{*}\,|\,x_{*},\boldsymbol{y},\boldsymbol{x},\boldsymbol{m})$
for new data $y_{*}$ given covariates $x_{*}$ and all training data. A
detailed explanation together with the necessary derivations to approximate
this intractable density is given in Appendix D.

##  4 Experiments

We demonstrate the efficiency of our proposal on various tasks, such as
missing value imputation, long-term prediction, for synthetic as well as real-
world healthcare datasets. We compare against a variety of models: GPPVAE
(Casale et al., 2018) serves as a general GP-prior representative, L-VAE
(Ramchandran et al., 2021) as a variant specifically designed for longitudinal
type of data, GRUI-GAN (Che et al., 2018) is a GAN-based model capable of
modelling non-random missingness, and mean imputation/prediction is a common
simple baseline. As GRUI-GAN is not designed for generative purposes, we only
provide imputation results for this method. For each method we evaluate mean-
squared error (MSE) and report mean performance as well as standard deviation
computed over five runs. The lowest mean in each experiment is marked bold in
the corresponding table. See Appendix E for further experimental details
(e.g., hyperparameters, kernel structures) that are not specified in the main
text and Appendix F for neural network architectures.

###  4.1 Regularly sampled Health MNIST

Table 2: Imputation MSEs on the regularly sampled health MNIST dataset. The
percentage represents maximum probability of pixel being missing.

Method 50% 75% 90% mean imputation
$0.266\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.314\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.373\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$ GPPVAE
$0.248\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.004$
$0.291\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.004$
$0.379\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.015$ GRUI-GAN
$0.224\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.037$
$0.218\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.012$
$0.269\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.028$ L-VAE
$\mathbf{0.124\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.009}$
$0.283\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.062$
$0.373\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$ LLSM (ours)
$\mathbf{0.124\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.008}$
$\mathbf{0.144\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.009}$
$\mathbf{0.174\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.016}$

Table 3: Future data prediction MSEs on the regularly sampled health MNIST
dataset. The percentage represents maximum probability of pixel being missing.

Method 50% 75% 90% mean prediction
$0.040\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.042\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.048\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$ GPPVAE
$0.041\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.041\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.048\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$ L-VAE
$\mathbf{0.021\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002}$
$0.038\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.008$
$0.047\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$ LLSM (ours)
$0.023\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002$
$\mathbf{0.024\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001}$
$\mathbf{0.026\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002}$

Figure 4: Future predictions of data (top) and missingness mask (bottom) on
the regularly sampled health MNIST dataset. Although predictions of mask look
almost identical, the prediction of data cannot be captured by L-VAE, whereas
LLSM does it very accurately.

To simulate a high-dimensional longitudinal dataset with structured
missingness, we used a modified version of the MNIST dataset (LeCun et al.,
1998) called Health MNIST (Krishnan et al., 2015). We chose two digits, ‘3’
and ‘6’ to represent two biological genders. We simulated a shared time-
related effect by shifting all digit instances towards the right corner over
time. In our experiments, half of the instances remain healthy and half get a
disease. To demonstrate changes in the laboratory measurements of the diseased
individuals, we rotated digits with the amount of rotation depending on the
time from disease diagnosis. Each sample has in total five covariates: _time_
, _id_ , _diseasePresence_ , _diseaseTime_ , and _gender_ , where _id_ stands
for the identifier of a specific instance. The timestamps of observations are
regularly sampled and hence only LLSM was evaluated.

To model MNAR, the probability of each pixel being missing depends on the
color intensity of that pixel: the higher the intensity, the higher is the
probability of the pixel being unobserved. For MAR, we applied a square-shaped
missingness mask to the images. If a patient is healthy, no box is applied.
For diseased patients, a mask is only applied upon the onset of the disease.
Afterwards, the mask starts to increase in size linearly as time progresses.
See Figure 4 for an illustration. In this case, missingness depends on both
unobserved signal as well as on covariate information.

The training set consists of $P=900$ unique instances, each having $n_{p}=20$
time points. The test set contains $100$ unique instances, with $15$ last
observations for each instance. When performing future prediction, the model
conditions on all training data as well as first five elements of each
instance, that are kept separately.

Table 2 shows that our method outperforms all baselines on the task of
missingness imputation. The same holds for future prediction of data from
covariates (Table 3), with the exception for the simplest missingness scenario
when L-VAE is slightly better. We also performed future prediction of
missingness mask by the same approach. Because none of the baseline methods is
able to model the mask $m$ explicitly, we separately modelled missingness by
training an L-VAE with a Bernoulli likelihood for $m$. The results are shown
in Table 4 and show that the mask prediction is almost identical except for
the case with the highest missingness when LLSM is slightly better. In Figure
4, we demonstrate visually the benefits of our model for future prediction for
the case of 75$\%$ maximum probability of missingness.

Table 4: Future missingness prediction MSEs on the regularly sampled health
MNIST dataset. The percentage represents maximum probability of pixel being
missing.

Method 50% 75% 90% L-VAE
$\mathbf{0.032\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000}$
$\mathbf{0.038\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002}$
$0.040\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002$ LLSM (ours)
$\mathbf{0.031\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002}$
$\mathbf{0.038\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.003}$
$\mathbf{0.038\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002}$

###  4.2 Irregularly sampled Health MNIST

Table 5: Imputation MSEs on the irregularly sampled health MNIST dataset. The
percentage represents maximum probability of pixel being missing.

Method 50% 75% 90% mean imputation
$0.259\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.335\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.380\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$ GPPVAE
$0.239\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002$
$0.319\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$
$0.396\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.004$ GRUI-GAN
$0.165\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.016$
$0.203\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.021$
$0.277\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.035$ L-VAE
$0.171\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.045$
$0.264\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.059$
$0.379\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$ LLSM (ours)
$0.130\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.007$
$0.163\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.011$
$\mathbf{0.191\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.012}$ LLPPSM (ours)
$\mathbf{0.128\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.005}$
$\mathbf{0.162\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.007}$
$0.207\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.028$

Figure 5: Future predictions of data (top) and missingness mask (bottom) on
irregularly sampled health MNIST dataset.

We modified the previous setup such that timepoints come from a random process
and similarity in the observations depends on covariates as well as the
underlying rate of the TPP that was discussed in Section 3.3. We implemented
it in a following way: for healthy patients, the timestamps come from a
homogeneous Poisson process with intensity $\lambda=0.1$ and the digit is not
modified, whereas for diseased patients, timestamps are generated according to
a Hawkes process, with baseline intensity $\mu=0.5$ and the digit rotation
depends on the intensity of process at the moment: the higher the intensity,
the stronger the rotation. For healthy patients we modelled MNAR as in the
previous setup, while for diseased patients we also applied a square-shaped
mask with its size being proportional to the intensity of the process at that
moment. Table 5 shows that both of our models improve upon the baselines in
all imputation scenarios. Moreover, Tables 6 and 7 show that LLPPSM
outperforms LLSM in both future data and mask prediction tasks. The Figure 5
depicts that although LLSM is capable of capturing the general form of an
image, it cannot model the rotation properly due to the limited kernel
component related to time whereas LLPPSM does it well. The same holds for
prediction of mask. The inferred mean intensity function of the point process
for one individual can be found in Appendix G.

###  4.3 Physionet data

We evaluated our model on healthcare data from the 2012 Physionet Challenge
(Silva et al., 2012). The dataset contains around 12,000 patients monitored on
the intensive care unit (ICU) for 48 hours. We modelled measurements of 37
different attributes, such as glucose level, heart rate, body temperature,
etc. The dataset is extremely sparse, with about 85$\%$ missing values. We
cannot directly measure imputation performance due to the lack of ground truth
data for missing values, hence, to test the learned representations of our
models, we used this dataset only for future prediction. We predicted values
of laboratory attributes given the knowledge of the first ten measurements for
a patient in the test set. As auxiliary covariates, we employed the following
variables: _time of the measurement_ , _id_ , _type of ICU_ , _gender_ , and
_in-hospital mortality_. More information regarding Physionet data can be
found in Section E.2. We present the results in Table 8. LLSM performs best,
with LLPPSM performing worse. This can be explained by the fact that many
observations are taken regularly each hour, making the temporal process be
pseudo-stochastic, hence modelling it just brings additional redundant
information to the model.

Table 6: Future data prediction MSEs on the irregularly sampled health MNIST
dataset. The percentage represents maximum probability of pixel being missing.

Method 50% 75% 90% mean prediction
$0.039\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.044\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.047\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$ GPPVAE
$0.040\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.044\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$
$0.049\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$ L-VAE
$0.029\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.005$
$0.036\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.006$
$0.047\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$ LLSM (ours)
$0.030\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$
$0.031\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$
$0.031\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$ LLPPSM (ours)
$\mathbf{0.025\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002}$
$\mathbf{0.027\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001}$
$\mathbf{0.030\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.003}$

Table 7: Future missingness prediction MSEs on the irregularly sampled health
MNIST dataset. The percentage represents maximum probability of pixel being
missing.

Method 50% 75% 90% L-VAE
$0.063\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$
$0.067\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$
$0.068\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$ LLSM (ours)
$0.063\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.002$
$0.067\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$
$0.066\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$ LLPPSM (ours)
$\mathbf{0.054\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.006}$
$\mathbf{0.056\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.005}$
$\mathbf{0.058\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.005}$

Table 8: Future prediction on Physionet dataset.

Method MSE mean prediction
$0.785\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.000$ GPPVAE
$0.786\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.001$ L-VAE
$0.720\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.008$ LLSM (ours)
$\mathbf{0.713\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.006}$ LLPPSM (ours)
$0.745\color[rgb]{0.2,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{0.2,0.2,0.2}\pgfsys@color@gray@stroke{0.2}\pgfsys@color<EMAIL_ADDRESS>0.009$

##  5 Conclusions.

In this work, we introduced a novel probabilistic framework for multivariate
data with missing values. First, we developed a deep latent variable model,
LLSM, that models structured missingness via separate set of latent variables.
Second, we extended this model by utilizing temporal point process to account
for stochastic nature of timepoints, LLPPSM. Our methods are specifically
designed for longitudinal type of data by leveraging GPs to define priors for
latent variables. We demonstrated excellent performance of both models on
different representation learning tasks and expect them to become useful tools
in the analysis of high-dimensional temporal and longitudinal data.

## Acknowledgements

We acknowledge the computational resources provided by the Aalto Science-IT
project. This work was supported by the Research Council of Finland (decision
number: 359135).

## Impact statement

This paper presents a work whose goal is to give a new idea of how high-
dimensional data with missing values can be analysed in the longitudinal
scenario. We do not see any potential harmful societal consequences of our
work.

## References

  * Casale et al. (2018) Casale, F. P., Dalca, A. V., Saglietti, L., Listgarten, J., and Fusi, N.  Gaussian process prior variational autoencoders.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Che et al. (2018) Che, Z., Purushotham, S., Cho, K., Sontag, D., and Liu, Y.  Recurrent neural networks for multivariate time series with missing value.  Scientific reports, 8(1):6085, 2018.
  * Collier et al. (2020) Collier, M., Nazabal, A., and Williams, C. K. I.  VAEs in the presence of missing data.  In _the First ICML Workshop on The Art of Learning with Missing Values_ , 2020.
  * Cox & Isham (1980) Cox, D. R. and Isham, V.  _Point processes_ , volume 12.  CRC Press, 1980.
  * Du et al. (2016) Du, N., Dai, H., Trivedi, R., Upadhyay, U., Gomez-Rodriguez, M., and Song, L.  Recurrent marked temporal point processes: Embedding event history to vector.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2016.
  * Fortuin et al. (2020) Fortuin, V., Baranchuk, D., Raetsch, G., and Mandt, S.  GP-VAE: Deep probabilistic time series imputation.  In _Proceedings of the 32rd International Conference on Artificial Intelligence and Statistics_ , 2020.
  * Ghalebikesabi et al. (2021) Ghalebikesabi, S., Cornish, R., Holmes, C., and Kelly, L.  Deep generative missingness pattern-set mixture models.  In _Proceedings of the 24th International Conference on Artificial Intelligence and Statistics_ , 2021.
  * Hawkes (1971) Hawkes, A. G.  Spectra of some self-exciting and mutually exciting point processes.  _Biometrika_ , 58(1):83–90, 1971.
  * Hawkes (2018) Hawkes, A. G.  Hawkes processes and their applications to finance: A review.  _Quant. Financ._ , 18(2):193–198, 2018.
  * Hensman et al. (2013) Hensman, J., Fusi, N., and Lawrence, N. D.  Gaussian processes for big data.  In _Proceedings of the 29th Conference on Uncertainty in Artificial Intelligence_ , 2013.
  * Ipsen et al. (2021) Ipsen, N. B., Mattei, P.-A., and Frellsen, J.  not-MIWAE: deep generative modelling with missing not at random data.  In _9th International Conference on Learning Representations_ , 2021.
  * John & Hensman (2018) John, S. and Hensman, J.  Large-scale cox process inference using variational fourier features.  In _Proceedings of the 35th International Conference on Machine Learning_ , 2018.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In _3rd International Conference on Learning Representations_ , 2015.
  * Kingma & Welling (2014) Kingma, D. P. and Welling, M.  Auto-encoding variational bayes.  In _2nd International Conference on Learning Representations_ , 2014.
  * Krishnan et al. (2015) Krishnan, R. G., Shalit, U., and Sontag, D.  Deep Kalman filters.  In _arXiv preprint arXiv:1511.05121_ , 2015.
  * LeCun et al. (1998) LeCun, Y., Cortes, C., and Burges, C.  The mnist database of handwritten digits.  1998\.
  * Liu & Hauskrecht (2019) Liu, S. and Hauskrecht, M.  Nonparametric regressive point processes based on conditional Gaussian processes.  In _Advances in Neural Information Processing Systems_ , 2019.
  * Liu (2015) Liu, X.  _Methods and applications of longitudinal data analysis_.  Elsevier, New York, 2015.
  * Lloyd et al. (2015) Lloyd, C., Gunter, T., Osborne, M., and Roberts, S.  Variational inference for Gaussian process modulated Poisson processes.  In _Proceedings of the 32nd International Conference on Machine Learning_ , 2015.
  * Luo et al. (2018) Luo, Y., Cai, X., ZHANG, Y., Xu, J., and xiaojie, Y.  Multivariate time series imputation with generative adversarial networks.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Mattei & Frellsen (2019) Mattei, P.-A. and Frellsen, J.  MIWAE: Deep generative modelling and imputation of incomplete data sets.  In _Proceedings of the 36th International Conference on Machine Learning_ , 2019.
  * Mitra et al. (2023) Mitra, R., McGough, S. F., Chakraborti, T., Holmes, C., Copping, R., Hagenbuch, N., Biedermann, S., Noonan, J., Lehmann, B., Shenvi, A., Doan, X. V., Leslie, D., Bianconi, G., Sanchez-Garcia, R., Davies, A., Mackintosh, M., Andrinopoulou, E.-R., Basiri, A., Harbron, C., and MacArthur, B. D.  Learning from data with structured missingness.  _Nat. Mach. Intell._ , 5:13–23, 2023.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems_ , 2019.
  * Puri & Tuan (1986) Puri, M. L. and Tuan, P. D.  Maximum likelihood estimation for stationary point processes.  In _Proceedings of the National Academy of Sciences_ , 1986.
  * Quiñonero Candela & Rasmussen (2005) Quiñonero Candela, J. and Rasmussen, C. E.  A unifying view of sparse approximate gaussian process regression.  In _Journal of Machine Learning Research_ , 2005.
  * Ramchandran et al. (2021) Ramchandran, S., Tikhonov, G., Kujanpää, K., Koskinen, M., and Lähdesmäki, H.  Longitudinal variational autoencoder.  In _Proceedings of the 24th International Conference on Artificial Intelligence and Statistics_ , 2021.
  * Rezende et al. (2014) Rezende, D. J., Mohamed, S., and Wierstra, D.  Stochastic backpropagation and approximate inference in deep generative models.  In _Proceedings of the 31st International Conference on Machine Learning_ , 2014.
  * Rubin (1976) Rubin, D. B.  Inference and missing data.  _Biometrika_ , 63:581–592, 1976.
  * Silva et al. (2012) Silva, I., Moody, G., Scott, D. J., Celi, L. A., and Mark, R. G.  Predicting in-hospital mortality of icu patients: The physionet/computing in cardiology challenge 2012.  In _2012 Computing in Cardiology.IEEE_ , 2012.
  * Titsias (2009) Titsias, M.  Variational learning of inducing variables in sparse gaussian processes.  In _Proceedings of the 12th International Conference on Artificial Intelligence and Statistics_ , 2009.

Appendix

##  Appendix A Deriving the ELBO for LLSM

By introducing a variational distribution
$q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})$,
the marginal likelihood can be expanded as

$\displaystyle\log
p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{x})$
$\displaystyle=\iint{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})\log{\frac{p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{x})}{p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}}d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}}$
$\displaystyle=\iint{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})\log{\frac{p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{m},\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x})}{p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}}d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}}$
$\displaystyle=\iint{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})\log{\frac{p_{\psi}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})p_{\theta}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x})}{p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}}d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}}$
$\displaystyle=\iint{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})\log{p_{\psi}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}}$
$\displaystyle\qquad+\iint{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})\log\frac{p_{\theta}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x})}{p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}}$
$\displaystyle={\operatorname*{\mathbb{E}}}_{q}\left[\log{p_{\psi}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}\right]$
$\displaystyle\qquad+\iint{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})\log\frac{p_{\theta}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x})q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}{p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}}$
$\displaystyle={\operatorname*{\mathbb{E}}}_{q}\left[\log{p_{\psi}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}\right]$
$\displaystyle\qquad+\iint{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})\log\frac{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}{p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}}$
$\displaystyle\qquad-\iint{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})\log\frac{q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})}{p_{\theta}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x})}d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}}$
$\displaystyle={\operatorname*{\mathbb{E}}}_{q}\left[\log{p_{\psi}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}\right]$
$\displaystyle\qquad+\text{KL}(q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})||p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m}))$
$\displaystyle\qquad-\text{KL}(q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})||p_{\theta}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x})).$

Hence,

$\displaystyle\log
p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{x})-\text{KL}(q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})||p_{\omega}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m}))=$
$\displaystyle\qquad{\operatorname*{\mathbb{E}}}_{q}\left[\log{p_{\psi}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}\right]-\text{KL}(q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})||p_{\theta}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x})).$

As the KL divergence term is positive, we get

$\displaystyle\log
p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{x})$
$\displaystyle\geq{\operatorname*{\mathbb{E}}}_{q}\left[\log{p_{\psi}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}\right]-\text{KL}(q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})||p_{\theta}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x}))$
$\displaystyle=\mathcal{L}(\phi,\psi,\theta;\boldsymbol{y}^{\text{o}},\boldsymbol{m}).$

By assuming the following factorizations:

$\displaystyle
p_{\psi}(\boldsymbol{y}^{\text{o}},\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})$
$\displaystyle=p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})$
$\displaystyle
q_{\phi}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{y}^{\text{o}},\boldsymbol{m})$
$\displaystyle=q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})$
$\displaystyle
p_{\theta}(\boldsymbol{z}^{y},\boldsymbol{z}^{m}|\boldsymbol{x})$
$\displaystyle=p_{\theta_{y}}(\boldsymbol{z}^{y}|\boldsymbol{x})p_{\theta_{m}}(\boldsymbol{z}^{m}|\boldsymbol{x}),$

the ELBO simplifies to

$\displaystyle\mathcal{L}(\phi,\psi,\theta;\boldsymbol{y}^{\text{o}},\boldsymbol{m})$
$\displaystyle={\operatorname*{\mathbb{E}}}_{q_{y}}\left[\log{p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})}\right]+{\operatorname*{\mathbb{E}}}_{q_{m}}\left[\log{p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}\right]$
$\displaystyle\qquad-\text{KL}(q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})||p_{\theta_{y}}(\boldsymbol{z}^{y}|\boldsymbol{x}))-\text{KL}(q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})||p_{\theta_{m}}(\boldsymbol{z}^{m}|\boldsymbol{x})),$

where $\psi_{y}$, $\psi_{m}$, $\phi_{y}$, $\phi_{m}$ are parameterized by
neural networks and $q_{y}$, $q_{m}$ are corresponding variational
distributions of $\boldsymbol{z}^{y}$ and $\boldsymbol{z}^{m}$.

By approximating the KL terms with the corresponding upper bounds
$\hat{{\text{KL}}}$ that are necessary for computational scalability (see
Appendix B), the final ELBO is given as

$\displaystyle\mathcal{L}(\phi,\psi,\theta;\boldsymbol{y}^{\text{o}},\boldsymbol{m})$
$\displaystyle\geq{\operatorname*{\mathbb{E}}}_{q_{y}}\left[\log{p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})}\right]+{\operatorname*{\mathbb{E}}}_{q_{m}}\left[\log{p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})}\right]$
$\displaystyle\qquad-\hat{{\text{KL}}}(q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})||p_{\theta_{y}}(\boldsymbol{z}^{y}|\boldsymbol{x}))-\hat{\text{KL}}(q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})||p_{\theta_{m}}(\boldsymbol{z}^{m}|\boldsymbol{x})).$

##  Appendix B Scalable KL Divergence Computation

Here we review the KL upper bound from (Ramchandran et al., 2021) that
implements a scalable KL divergence computation. Optimising the variational
objective requires us to evaluate

$\displaystyle\text{KL}(q_{\phi}(\boldsymbol{z}|\boldsymbol{y})||p_{\theta}(\boldsymbol{z}|\boldsymbol{x}))=\text{KL}({\mathcal{N}(\bar{\boldsymbol{\mu}},W)||\mathcal{N}(\boldsymbol{0},\Sigma))},$
(5)

where
$\bar{\boldsymbol{\mu}}=[\mu_{\phi}(y_{1}),\ldots,\mu_{\phi}(y_{N})]^{T},W=\mathrm{diag}(\sigma^{2}_{\phi}(y_{1}),\ldots,\sigma^{2}_{\phi}(y_{N}))$,
and
$\Sigma=\sum_{r=1}^{R}{K^{(r)}_{\boldsymbol{x}\boldsymbol{x}}+\sigma_{zl}^{2}I_{N}}$.
The exact computation has $\mathcal{O}(N^{3})$ complexity, making it
impractical for large datasets. Therefore, instead of computing it, we will
use an upper bound to the KL that comes from the fact that any lower bound for
the prior GP marginal log-likelihood induces an upper bound to the KL
divergence. Titsias (2009) proposed the free-form variational lower bound for
a GP marginal log-likelihood
$\log\mathcal{N}(\boldsymbol{z}\big{|}\boldsymbol{0},\Sigma)$ by assuming a
set of $M$ inducing locations $\boldsymbol{s}=[s_{1},\ldots,s_{M}]$ in
$\mathcal{X}$, with the corresponding inducing function values
$\boldsymbol{u}=[f(s_{1}),\ldots,f(s_{M})]^{T}=[u_{1},\ldots,u_{M}]^{T}$, such
that

$\displaystyle p(\boldsymbol{u})$
$\displaystyle=\mathcal{N}(\boldsymbol{u}\big{|}\boldsymbol{0},K_{\boldsymbol{s}\boldsymbol{s}})$
$\displaystyle p(\boldsymbol{f}\big{|}\boldsymbol{u})$
$\displaystyle=\mathcal{N}(\boldsymbol{f}\big{|}K_{\boldsymbol{x}\boldsymbol{s}}K^{-1}_{\boldsymbol{s}\boldsymbol{s}}\boldsymbol{u},\tilde{K}),$
$\displaystyle\tilde{K}$
$\displaystyle=K_{\boldsymbol{x}\boldsymbol{x}}-K_{\boldsymbol{x}\boldsymbol{s}}K^{-1}_{\boldsymbol{s}\boldsymbol{s}}K_{\boldsymbol{s}\boldsymbol{x}}$
$\displaystyle p(\boldsymbol{z}\big{|}\boldsymbol{f})$
$\displaystyle=\mathcal{N}(\boldsymbol{z}\big{|}\boldsymbol{f},\sigma^{2}_{z}I_{N}),$

and the corresponding lower bound is
$\mathcal{L}(\boldsymbol{z};\Sigma)=\mathcal{N}(\boldsymbol{z}\big{|}\boldsymbol{0},K_{\boldsymbol{x}\boldsymbol{s}}K^{-1}_{\boldsymbol{s}\boldsymbol{s}}K_{\boldsymbol{s}\boldsymbol{x}}+\sigma^{2}_{z}I_{N})-\frac{1}{2\sigma^{2}_{z}}\text{tr}(\tilde{K})$,
where $\text{tr}(\cdot)$ is the matrix trace. Titsias bound is known to be
tight when $M$ is large enough and the covariance function is smooth enough.
Longitudinal data, however, always contain categorical covariates
corresponding to instance _ids_ , making the covariance function necessarily
non-continuous.

To still get a tighter bound, we separate the additive component that
corresponds to the instance _id_ from the other covariates, so that covariance
function has the following form
$\Sigma=K^{(A)}_{\boldsymbol{x}\boldsymbol{x}}+\hat{\Sigma}$, where
$\hat{\Sigma}=\mathrm{diag}(\hat{\Sigma}_{1},\ldots,\hat{\Sigma}_{P})$,
$\hat{\Sigma}_{p}=K^{(R)}_{\boldsymbol{x}_{p}\boldsymbol{x}_{p}}+\sigma^{2}_{z}I_{n_{p}}$
contains all terms with instance-specific _id_ and
$K^{(A)}_{\boldsymbol{x}\boldsymbol{x}}=\sum_{r=1}^{R-1}{K^{(r)}_{\boldsymbol{x}\boldsymbol{x}}}$
contains the other components. Ramchandran et al. (2021) proposed the
following upper bound for KL

$\displaystyle\text{KL}\leq\frac{1}{2}\left(\text{tr}\left(\bar{\Sigma}^{-1}W\right)+\bar{\boldsymbol{\mu}}^{T}\bar{\Sigma}^{-1}\bar{\boldsymbol{\mu}}-N+\log|\bar{\Sigma}|-\log|W|+\sum_{p=1}^{P}{\text{tr}\left(\hat{\Sigma}^{-1}_{p}\tilde{K}^{(A)}_{\boldsymbol{x}_{p}\boldsymbol{x}_{p}}\right)}\right),$

where
$\bar{\Sigma}=K^{(A)}_{\boldsymbol{x}\boldsymbol{s}}K^{{(A)}^{-1}}_{\boldsymbol{s}\boldsymbol{s}}K^{(A)}_{\boldsymbol{s}\boldsymbol{x}}+\hat{\Sigma}$
and
$\tilde{K}^{(A)}_{\boldsymbol{x}_{p}\boldsymbol{x}_{p}}=K^{(A)}_{\boldsymbol{x}_{p}\boldsymbol{x}_{p}}-K^{(A)}_{\boldsymbol{x}_{p}\boldsymbol{s}}K^{{(A)}^{-1}}_{\boldsymbol{s}\boldsymbol{s}}K^{(A)}_{\boldsymbol{s}\boldsymbol{x}_{p}}$.
This bound has a computational complexity
$\mathcal{O}(\sum_{p=1}^{P}{{n_{p}}^{3}+NM^{2})}$ and Ramchandran et al.
(2021) proved it to be tighter than the corresponding bound by Titsias for
longitudinal datasets. Despite the reduced complexity, a problem still
remains. For every gradient descent step, the algorithm has to iterate through
the entire dataset, requiring a substantial allocation of memory and
computational time. This issue can be solved using a similar technique as the
one presented by Hensman et al. (2013), with adaptation to the properties of a
longitudinal kernel. We will only present the final bound and refer the reader
to Ramchandran et al. (2021) for a detailed derivation. Defining $I_{p_{i}}$
to be the index of the $i$th sample for the $p$th patient and
$\bar{\boldsymbol{\mu}}_{p}=\left[\bar{\boldsymbol{\mu}}_{I_{p_{i}}},\ldots,\bar{\boldsymbol{\mu}}_{I_{n_{p}}}\right]^{T}$
to be a a sub-vector of $\bar{\boldsymbol{\mu}}$ that is related to the $p$th
patient, the unbiased estimate of the KL divergence upper bound, computed from
the batch with instances $\mathrm{P}_{\text{batch}}\subset{\\{1,\ldots,P\\}}$,
is defined as

$\displaystyle\hat{\text{KL}}$
$\displaystyle=\frac{1}{2}\frac{P}{|\mathrm{P}_{\text{batch}}|}\sum_{p\in\mathrm{P}}\left(\left(K^{(A)}_{\boldsymbol{x}_{p}\boldsymbol{s}}K^{{(A)}^{-1}}_{\boldsymbol{s}\boldsymbol{s}}\boldsymbol{m}_{H}-\bar{\boldsymbol{\mu}}_{p}\right)^{T}{\hat{\Sigma}^{-1}}_{p}\left(K^{(A)}_{\boldsymbol{x}_{p}\boldsymbol{s}}K^{{(A)}^{-1}}_{\boldsymbol{s}\boldsymbol{s}}\boldsymbol{m}_{H}-\bar{\boldsymbol{\mu}}_{p}\right)+\sum_{i=1}^{n_{p}}({\hat{\Sigma}}^{-1}_{p})_{ii}\sigma^{2}_{\phi}(\boldsymbol{y}_{I_{p_{i}}})\right.$
$\displaystyle\qquad+\left.\quad\log|\hat{\Sigma}_{p}|+\text{tr}\left(\hat{\Sigma}^{-1}_{p}{\tilde{K}^{(A)}}_{\boldsymbol{x}_{p}\boldsymbol{x}_{p}}\right)+\text{tr}\left(\left(K^{{(A)}^{-1}}_{\boldsymbol{s}\boldsymbol{s}}HK^{{(A)}^{-1}}_{\boldsymbol{s}\boldsymbol{s}}\right)\left(K^{(A)}_{\boldsymbol{s}\boldsymbol{x}_{p}}\hat{\Sigma}^{-1}_{p}K^{(A)}_{\boldsymbol{x}_{p}\boldsymbol{s}}\right)\right)\right.\left.-\sum_{i=1}^{n_{p}}\log\sigma^{2}_{\phi}(\boldsymbol{y}_{I_{p_{i}}})\right)$
$\displaystyle\qquad-\frac{N}{2}+\text{KL}\left[\mathcal{N}(\boldsymbol{m}_{H},H)||\mathcal{N}(\boldsymbol{0},K^{(A)}_{\boldsymbol{s}\boldsymbol{s}})\right],$

where $\boldsymbol{m}_{H}$ and $H$ are variational parameters computed via
natural gradients.

##  Appendix C Various LLPPSM specifications

We define three sets of latent variables: $z^{y},z^{m}$ and $z^{\lambda}$.
Observations $\boldsymbol{y}$, masks $\boldsymbol{m}$ and timestamps
$\boldsymbol{t}$ are modelled as random variables. The complete joint
probability is given as:

$\displaystyle
p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{z}^{y},\boldsymbol{m},\boldsymbol{z}^{m},\boldsymbol{t},z^{\lambda}|\boldsymbol{x}^{\text{s}},\boldsymbol{v}_{D})=$
$\displaystyle
p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})p_{\theta_{y}}(\boldsymbol{z}^{y}|\boldsymbol{x},\lambda(\boldsymbol{t}))p_{\theta_{m}}(\boldsymbol{z}^{m}|\boldsymbol{x},\lambda(\boldsymbol{t}))$
$\displaystyle\qquad
p(\boldsymbol{t}|\lambda)p_{\theta_{\lambda}}(z^{\lambda}|\boldsymbol{v}_{D}),$

where $\boldsymbol{v}_{D}$ corresponds to elapsed times between
$\boldsymbol{t}$ and $D$ previous timepoints that occurred before
$\boldsymbol{t}$, denoted as $\boldsymbol{t}_{D}$, $\boldsymbol{x}^{\text{s}}$
refers to static covariates and $\boldsymbol{x}$ consists of $\boldsymbol{t}$
and $\boldsymbol{x}^{\text{s}}$. To compute the marginal likelihood, we have
to marginalize over the latent variables as

$\displaystyle
p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{m},\boldsymbol{t}|\boldsymbol{x}^{\text{s}},\boldsymbol{v}_{D})$
$\displaystyle=\iiint{p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{z}^{y},\boldsymbol{m},\boldsymbol{z}^{m},\boldsymbol{t},z^{\lambda}|\boldsymbol{x}^{\text{s}},\boldsymbol{v}_{D})d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}dz^{\lambda}}$
$\displaystyle=\iiint{p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})p_{\theta_{y}}(\boldsymbol{z}^{y}|\boldsymbol{x},\lambda(\boldsymbol{t}))p_{\theta_{m}}(\boldsymbol{z}^{m}|\boldsymbol{x},\lambda(\boldsymbol{t}))}$
$\displaystyle\qquad\cdot
p(\boldsymbol{t}|\lambda)p_{\theta_{\lambda}}(z^{\lambda}|\boldsymbol{v}_{D})d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}dz^{\lambda}.$

Factorizing the joint likelihood w.r.t. the observation given the latent
variables, we get

$\displaystyle
p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{m},\boldsymbol{t}|\boldsymbol{x}^{\text{s}},\boldsymbol{v}_{D})=\iiint$
$\displaystyle\prod_{i=1}^{N}{p_{\psi_{y}}(y^{\text{o}}_{i}|z_{i}^{y},z_{i}^{m},m_{i})p_{\psi_{m}}(m_{i}|z_{i}^{y},z_{i}^{m})}p_{\theta_{y}}(\boldsymbol{z}^{y}|\boldsymbol{x},\lambda(\boldsymbol{t}))p_{\theta_{m}}(\boldsymbol{z}^{m}|\boldsymbol{x},\lambda(\boldsymbol{t}))$
$\displaystyle\qquad
p(\boldsymbol{t}|\lambda)p_{\theta_{\lambda}}(z^{\lambda}|\boldsymbol{v}_{D})d\boldsymbol{z}^{y}d\boldsymbol{z}^{m}dz^{\lambda}.$

We assume the following Gaussian process priors

$\displaystyle z^{\lambda}(t)$
$\displaystyle\sim{\mathcal{G}}P\left(\big{.}0,k(v_{D},{v_{D}}^{\prime}|\theta_{\lambda})\right)$
$\displaystyle z^{y}(x,\lambda)$
$\displaystyle\sim{\mathcal{G}}P\left(\big{.}0,k((x,\lambda(t)),(x^{\prime},\lambda(t^{\prime}))|\theta_{y})\right)$
$\displaystyle z^{m}(x,\lambda)$
$\displaystyle\sim{\mathcal{G}}P\left(\big{.}0,k((x,\lambda(t)),(x^{\prime},\lambda(t^{\prime}))|\theta_{m})\right).$

The decoders of $\boldsymbol{y}$ and $\boldsymbol{m}$ are parameterized by
neural networks that predict the mean of the generative distributions

$\displaystyle p_{\psi_{y}}(y^{\text{o}}_{i}|z_{i}^{y},z_{i}^{m},m_{i})$
$\displaystyle=N(y_{i}\big{|}g_{\psi_{y}}(z_{i}^{y},z_{i}^{m}),\text{diag}(\sigma^{2}_{y_{1}},\ldots,\sigma^{2}_{y_{K}}))\odot
m_{i}$ $\displaystyle p_{\psi_{m}}(m_{i}\big{|}z_{i}^{y},z_{i}^{m})$
$\displaystyle=Ber(m_{i}|g_{\psi_{m}}(z_{i}^{y},z_{i}^{m})),$

where $\sigma^{2}_{y_{1}},\ldots,\sigma^{2}_{y_{K}}$ are jointly optimised via
gradient descent, and $\odot$ is the elementwise Hadamard product.

The likelihood of the temporal point process for instance $p$ is given as

$\displaystyle p(\boldsymbol{t}_{p}|\lambda)$
$\displaystyle=\exp\Big{(}-\int_{\mathcal{T}}{\lambda(t)dt}\Big{)}\prod_{i=1}^{n_{p}}{\lambda(t_{i})},$
$\displaystyle\text{where}\qquad\lambda(t)$
$\displaystyle=(z^{\lambda}(t)+\beta)^{2},$

and $\beta$ is a learnable offset.

The covariances of the GPs of $z^{y}$ and $z^{m}$ are parameterized by the
additive kernels discussed in Section 3. For $z^{\lambda}$, we instead rely on
the following additive structure:

$\displaystyle k(v_{D},{v_{D}}^{\prime}|\theta_{\lambda})$
$\displaystyle=\sum_{d=1}^{D}{\mathbbm{1}({{v_{d}}})\mathbbm{1}({{{v_{d}}^{\prime}}})\cdot\gamma_{d}\cdot\exp\left(-\frac{(v_{d}-{v_{d}}^{{}^{\prime}})^{2}}{2l_{d}^{2}}\right)},$
$\displaystyle\mathbbm{1}({{v_{d}}})$ $\displaystyle=\begin{cases}1,\text{if}\
{v_{d}}<\infty\\\ 0,\text{otherwise}\\\ \end{cases},$

where $v_{d}$ is the elapsed time between $t$ and $d$th timepoint, $t_{d}$,
that happened before $t$, and infinite, if there is no available information
about the past event. Hence, the above kernel depends on past events and on
the current time point.

As the posterior inference is not tractable, we rely on variational inference
and choose the following approximation to the posterior

$q(\boldsymbol{z}^{y},\boldsymbol{z}^{m},z^{\lambda},\boldsymbol{u}\big{|}\boldsymbol{y}^{\text{o}},\boldsymbol{m})=q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})\cdot
q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})\cdot
p(z^{\lambda}|\boldsymbol{u})\cdot q(\boldsymbol{u}),$

where we also employed inducing points of $z^{\lambda}$, denoted by
$\boldsymbol{u}$ with $q(\boldsymbol{u})=N(m_{\lambda},S)$. Note that these
inducing points are different from those of the scalable KL bound in Appendix
B. To shorten the notation we will denote this variational posterior as $q$ in
the following derivations. The ELBO is given as

$\displaystyle\mathcal{L}$
$\displaystyle={\operatorname*{\mathbb{E}}}_{q}\left[\log\frac{p_{\omega}(\boldsymbol{y}^{\text{o}},\boldsymbol{m},\boldsymbol{t},\boldsymbol{z}^{y},\boldsymbol{z}^{m},z^{\lambda},\boldsymbol{u}|\boldsymbol{x}^{\text{s}},\boldsymbol{v}_{D})}{q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})\cdot
q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})\cdot
p(z^{\lambda}|\boldsymbol{u})\cdot q(\boldsymbol{u})}\right]$
$\displaystyle={\operatorname*{\mathbb{E}}}_{q}\left[\log\left(p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})p(\boldsymbol{t}|\lambda)\right)\right]-{\operatorname*{\mathbb{E}}}_{q}\left[\log(\frac{q(\boldsymbol{u})}{p(\boldsymbol{u})})\right]$
$\displaystyle\qquad-{\operatorname*{\mathbb{E}}}_{q}\left[\log(\frac{q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})}{p(\boldsymbol{z}^{y}|\boldsymbol{x},\lambda(\boldsymbol{t}))})\right]-{\operatorname*{\mathbb{E}}}_{q}\left[\log(\frac{q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})}{p(\boldsymbol{z}^{m}|\boldsymbol{x},\lambda(\boldsymbol{t}))})\right]$
$\displaystyle={\operatorname*{\mathbb{E}}}_{q}\left[\log\left(p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})p(\boldsymbol{t}|\lambda)\right)\right]-\text{KL}(q(\boldsymbol{u})||p(\boldsymbol{u}))$
$\displaystyle\qquad-{\operatorname*{\mathbb{E}}}_{p(z^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})}\left[\text{KL}(q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})||p(\boldsymbol{z}^{y}|\boldsymbol{x},\lambda(\boldsymbol{t}))\right]-{\operatorname*{\mathbb{E}}}_{p(z^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})}\left[\text{KL}(q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})||p(\boldsymbol{z}^{m}|\boldsymbol{x},\lambda(\boldsymbol{t}))\right]$
$\displaystyle\geq{\operatorname*{\mathbb{E}}}_{q}\left[\log\left(p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})p(\boldsymbol{t}|\lambda)\right)\right]-\text{KL}(q(\boldsymbol{u})||p(\boldsymbol{u}))$
$\displaystyle\qquad-{\operatorname*{\mathbb{E}}}_{p(z^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})}\left[\hat{\text{KL}}(q_{\phi_{y}}(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})||p(\boldsymbol{z}^{y}|\boldsymbol{x},\lambda(\boldsymbol{t}))\right]-{\operatorname*{\mathbb{E}}}_{p(z^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})}\left[\hat{\text{KL}}(q_{\phi_{m}}(\boldsymbol{z}^{m}|\boldsymbol{m})||p(\boldsymbol{z}^{m}|\boldsymbol{x},\lambda(\boldsymbol{t}))\right],$

where $\hat{\text{KL}}$ is the corresponding upper bound for KL divergence
from Appendix B and expectations involving upper bounds are estimated via
Monte Carlo sampling. The generative part can be further decomposed into

$\displaystyle{\operatorname*{\mathbb{E}}}_{q}\left[\log\left(p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})p(\boldsymbol{t}|\lambda)\right)\right]$
$\displaystyle\qquad={\operatorname*{\mathbb{E}}}_{q_{\phi_{y}}{(\boldsymbol{z}^{y}|\boldsymbol{y}^{\text{o}})}}{\left[\log\left(p_{\psi_{y}}(\boldsymbol{y}^{\text{o}}|\boldsymbol{z}^{y},\boldsymbol{z}^{m},\boldsymbol{m})\right)\right]}+{\operatorname*{\mathbb{E}}}_{q_{\phi_{m}}{(\boldsymbol{z}^{m}|\boldsymbol{m})}}{\left[\log\left(p_{\psi_{m}}(\boldsymbol{m}|\boldsymbol{z}^{y},\boldsymbol{z}^{m})\right)\right]}+{\operatorname*{\mathbb{E}}}_{p(z^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})}\left[\log(p(\boldsymbol{t}|\lambda)\right].$

The expectations involving $\boldsymbol{z}^{y}$ and $\boldsymbol{z}^{m}$ are
computed by sampling from their variational distributions and utilizing the
reparameterization trick (Kingma & Welling, 2014), whereas expectation
involving likelihood of the point process can be evaluated in a closed from
due to the chosen squared link function as we show below.

To lighten the notation, we use
$\mathcal{L}_{T}:={\operatorname*{\mathbb{E}}}_{p(\boldsymbol{z}^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})}\left[\log(p(\boldsymbol{t}|\lambda)\right]$, drop subscript
$p$ denoting the individual and drop subscript $D$ from $v_{D}$. Also, by
$\boldsymbol{s}$ we denote inducing point locations of $\boldsymbol{u}$.
First, we integrate out $\boldsymbol{u}$:

$\displaystyle q(z^{\lambda})$
$\displaystyle=\int_{U}{p(z^{\lambda}|\boldsymbol{u})\cdot
q(\boldsymbol{u})d\boldsymbol{u}}=N(z^{\lambda}\big{|}\tilde{\mu},\tilde{\Sigma}),$
$\displaystyle\text{where}\quad\tilde{\mu}(v)$
$\displaystyle=k_{v\boldsymbol{s}}K_{\boldsymbol{s}\boldsymbol{s}}^{-1}m_{\lambda},$
$\displaystyle\tilde{\Sigma}(v,{v}^{\prime})$
$\displaystyle=K_{v{v}^{\prime}}-k_{v\boldsymbol{s}}K_{\boldsymbol{s}\boldsymbol{s}}^{-1}k_{\boldsymbol{s}{v}^{\prime}}+k_{v\boldsymbol{s}}K_{\boldsymbol{s}\boldsymbol{s}}^{-1}SK_{\boldsymbol{s}\boldsymbol{s}}^{-1}k_{\boldsymbol{s}{v}^{\prime}},$

and $U$ is a space of inducing values. We write $\mathcal{L}_{T}$ as

$\mathcal{L}_{T}=\sum_{n}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[\log\lambda(t_{n})\right]}-{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\Big{[}\int_{T}{\lambda(t)dt}\Big{]}=\sum_{n}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[\log\left({z^{\lambda}(t_{n})+\beta}\right)^{2}\right]}-\underbrace{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\Big{[}\int_{T}{\lambda(t)dt}\Big{]}}_{:=\mathcal{L}_{t}}$

where we sum over all timepoints of the individual. Via Lloyd et al. (2015),
we have

$\displaystyle{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[\log\left({z^{\lambda}(t_{n})+\beta}\right)^{2}\right]$
$\displaystyle=\int_{-\infty}^{\infty}{\log({(z^{\lambda}(t_{n})+\beta})^{2})N(z^{\lambda}(t_{n})\big{|}\tilde{\mu},\tilde{\sigma}^{2})dz^{\lambda}(t_{n})}$
$\displaystyle=-\tilde{G}\Big{(}-\frac{({\tilde{\mu}+\beta})^{2}}{2\tilde{\sigma}^{2}}\Big{)}+\log\Big{(}\frac{\tilde{\sigma}^{2}}{2}\Big{)}-C.$

where $\tilde{\sigma}^{2}$ is the diagonal element of $\tilde{\Sigma}$, $C$ is
the Euler-Mascheroni constant and $\tilde{G}$ is the confluent hyper-geometric
function.

Following the derivations of Liu & Hauskrecht (2019), we compute
$\mathcal{L}_{t}$ as:

$\displaystyle\mathcal{L}_{t}$
$\displaystyle={\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[\int_{T}{\lambda(t)dt}\right]={\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[\int_{T}{(z^{\lambda}(t)+\beta)^{2}dt}\right]$
$\displaystyle=\int_{T}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[(z^{\lambda}(t)+\beta)^{2}\right]dt}$
$\displaystyle=\int_{T}{\left({\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[(z^{\lambda}(t)^{2}\right]+2\beta{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]+\beta^{2}\right)dt}$
$\displaystyle=\int_{T}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]^{2}dt}+\int_{T}{\text{Var}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]dt}+2\beta\int_{T}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]dt}+\beta^{2}|T|$
$\displaystyle=\sum_{n}{\left[\int_{t_{n-1}}^{t_{n}}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]^{2}dt}+\int_{t_{n-1}}^{t_{n}}{\text{Var}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]dt}+2\beta\int_{t_{n-1}}^{t_{n}}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]dt}\right]+\beta^{2}|T|}.$

Each integral is computed as follows:

$\displaystyle\int_{\boldsymbol{t}_{n-1}}^{\boldsymbol{t}_{n}}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]^{2}dt}$
$\displaystyle=m_{\lambda}^{T}K_{\boldsymbol{s}\boldsymbol{s}}^{-1}\Psi_{n}K_{\boldsymbol{s}\boldsymbol{s}}^{-1}m_{\lambda}$
$\displaystyle\int_{t_{n-1}}^{t_{n}}{\text{Var}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]dt}$
$\displaystyle=\sum_{d=1}^{D}{\gamma_{d}\int_{t_{n-1}}^{t_{n}}{\mathbbm{1}}({v_{d}})dt}-\text{tr}(K_{\boldsymbol{s}\boldsymbol{s}}^{-1}\Psi_{n})+\text{tr}(K_{\boldsymbol{s}\boldsymbol{s}}^{-1}SK_{\boldsymbol{s}\boldsymbol{s}}^{-1}\Psi_{n})$
$\displaystyle\int_{t_{n-1}}^{t_{n}}{{\operatorname*{\mathbb{E}}}_{q(z^{\lambda})}\left[z^{\lambda}(t)\right]dt}$
$\displaystyle=\Phi_{n}(\boldsymbol{s})^{T}K_{\boldsymbol{s}\boldsymbol{s}}^{-1}m_{\lambda}$
$\displaystyle\Psi_{n}(\boldsymbol{s},\boldsymbol{s}^{\prime})$
$\displaystyle=\int_{t_{n-1}}^{t_{n}}{K(\boldsymbol{s},v(t))K(v(t),\boldsymbol{s}^{\prime})dt}$
$\displaystyle\Phi_{n}(\boldsymbol{s})$
$\displaystyle=\int_{t_{n-1}}^{t_{n}}{K(\boldsymbol{s},v(t))dt}.$

$\Phi$ and $\Psi$ each have closed form solutions which we obtain by
evaluating the integrals for the sum of SE kernels.

$\displaystyle\Phi_{n}(\boldsymbol{s})$
$\displaystyle=\sum_{d=1}^{D}{\mathbbm{1}(\boldsymbol{s}_{d})\mathbbm{1}(v_{d_{n}})\gamma_{d}\frac{\sqrt{\pi
l_{d}^{2}}}{\sqrt{2}}\left[\text{erf}\left(\frac{{v_{d}}_{n}-\boldsymbol{s}_{d}}{\sqrt{2l_{d}^{2}}}\right)-\text{erf}\left(\frac{{v_{d}}_{n-1}-\boldsymbol{s}_{d}}{\sqrt{2l_{d}^{2}}}\right)\right]}$
$\displaystyle\Psi_{n}(\boldsymbol{s},\boldsymbol{s}^{\prime})$
$\displaystyle=\sum_{i,j=1}^{D}\mathbbm{1}(\boldsymbol{s}_{i})\mathbbm{1}(\boldsymbol{s}_{j}^{\prime})\mathbbm{1}({v_{i}}_{n})\mathbbm{1}({v_{j}}_{n})\gamma_{i}\gamma_{j}\frac{\sqrt{\pi
l_{i}^{2}l_{j}^{2}}}{\sqrt{2\cdot(l_{i}^{2}+l_{j}^{2})}}\exp\left(-\frac{(\boldsymbol{s}_{i}+{t_{i}}_{n}-\boldsymbol{s}_{j}^{\prime}-{t_{j}}_{n})^{2}}{2(l_{i}^{2}+l_{j}^{2})}\right)$
$\displaystyle\qquad\left[\text{erf}\left(\frac{l_{i}^{2}({v_{j}}_{n}-\boldsymbol{s}_{j}^{\prime})+l_{j}^{2}({v_{i}}_{n}-\boldsymbol{s}_{i})}{\sqrt{2l_{i}^{2}l_{j}^{2}(l_{i}^{2}+l_{j}^{2})}}\right)-\text{erf}\left(\frac{l_{i}^{2}({v_{j}}_{n-1}-\boldsymbol{s}_{j}^{\prime})+l_{j}^{2}({v_{i}}_{n-1}-\boldsymbol{s}_{i})}{\sqrt{2l_{i}^{2}l_{j}^{2}(l_{i}^{2}+l_{j}^{2})}}\right)\right],$

where $l_{i}$, $\gamma_{i}$ are kernel hyperparameters.

##  Appendix D Predictive distribution

Given training samples $\boldsymbol{y}^{\text{o}}$, covariate information
$\boldsymbol{x}$ and masks $\boldsymbol{m}$, the predictive distribution for a
new observation $y_{*}$, given covariates $x_{*}$ is

$\displaystyle
p(y_{*}\,|\,x_{*},\boldsymbol{y}^{\text{o}},\boldsymbol{x},\boldsymbol{m})$
$\displaystyle=\int
p(y_{*}\,|\,z^{y}_{*},z^{m}_{*})p(z^{y}_{*},z^{m}_{*}\,|\,x_{*},\boldsymbol{y}^{\text{o}},\boldsymbol{x},\boldsymbol{m})\,dz^{y}_{*}\,dz^{m}_{*}$
$\displaystyle\approx\int\underbrace{p(y_{*}\,|\,z^{y}_{*},z^{m}_{*})}_{\text{decode
GP
predictions}}\underbrace{p(z^{y}_{*}\,|\,x_{*},\lambda(t_{*}),\boldsymbol{z}^{y},\boldsymbol{x},\lambda(\boldsymbol{t}))}_{\text{GP
posterior of
$z^{y}_{*}$}}\underbrace{p(z^{m}_{*}\,|\,x_{*},\lambda(t_{*}),\boldsymbol{z}^{m},\boldsymbol{x},\lambda(\boldsymbol{t}))}_{\text{GP
posterior of $z^{m}_{*}$}}$
$\displaystyle\quad\underbrace{q(z^{\lambda}_{*})}_{\text{variational
posterior of
$z^{\lambda}_{*}$}}\underbrace{q_{\phi_{y}}(\boldsymbol{z}^{y}\,|\,\boldsymbol{y}^{\text{o}})}_{\text{encode
data}}\underbrace{q_{\phi_{m}}(\boldsymbol{z}^{m}\,|\,\boldsymbol{m})}_{\text{encode
mask}}\,dz^{\lambda}_{*}\,dz^{y}_{*}\,dz^{m}_{*}\,d\boldsymbol{z}^{y}\,d\boldsymbol{z}^{m},$

where timestamps $\boldsymbol{t}$ and $t_{*}$ belong to $\boldsymbol{x}$ and
$x_{*}$, respectively, and
$p(z^{y}_{*}\,|\,x_{*},\lambda(t_{*}),\boldsymbol{z}^{y},\boldsymbol{x},\lambda(\boldsymbol{t}))$
and
$p(z^{m}_{*}\,|\,x_{*},\lambda(t_{*}),\boldsymbol{z}^{m},\boldsymbol{x},\lambda(\boldsymbol{t}))$
are GP posteriors such that

$\displaystyle
p(z_{*}\,|\,x_{*},\lambda(t_{*}),\boldsymbol{z},\boldsymbol{x},\lambda(\boldsymbol{t}))=\mathcal{N}(\tilde{\mu},\tilde{\sigma}^{2}),$
$\displaystyle\tilde{\mu}={k_{*}}^{T}(K_{\boldsymbol{x}\boldsymbol{x}}+\sigma_{z}^{2}I_{N})^{-1}\boldsymbol{z},$
$\displaystyle\tilde{\sigma}^{2}=k_{*}+\sigma_{z}^{2}-{k}^{T}(K_{\boldsymbol{x}\boldsymbol{x}}+\sigma_{z}^{2}I_{N})^{-1}k,$
$\displaystyle k_{*}=k((x_{*},\lambda(t_{*})),(x_{*},\lambda(t_{*}))|\theta),$
$\displaystyle
k=k((x_{*},\lambda(t_{*})),(\boldsymbol{x},\lambda(\boldsymbol{t}))|\theta).$

Above, we dropped the superscripts, meaning that the same formulae hold for
both $z^{y}$ and $z^{m}$ with respect to their kernel hyperparameters. The
same approach holds for deriving
$p(m_{*}\,|\,x_{*},\boldsymbol{y}^{\text{o}},\boldsymbol{x},\boldsymbol{m})$.
In order to sample from these predictive distributions, ancestral sampling can
be employed. Computing the predictive distributions scales cubically for this
case. To get the idea of how to implement a scalable predictive distribution
using low-rank inducing point approximation, we refer the reader to the
derivations by Ramchandran et al. (2021).

##  Appendix E Experimental setups

We employed identical kernel structures for the GPs of both $z^{y}$ and
$z^{m}$ across all of the datasets mentioned below. Nonetheless, it’s
important to note that, in general, these kernel structures could differ
between the two. We also selected sixty inducing points for each GP model for
all setups and chose the latent dimension to be 32.

###  E.1 HealthMNIST variants

For the _Health MNIST regularly sampled_ dataset we use the following
covariates: _time_ , _id_ , _diseasePresence_ , _diseaseTime_ and _gender_.
When running LLSM on the corresponding dataset, we relied on the following
additive kernel structure

$f_{\text{ca}}(\text{id})+f_{\text{se}}(\text{time})+f_{\text{ca}\times\text{se}}(\text{id}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{gender}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{diseasePresence}\times\text{diseaseTime}),$

where se denotes squared exponential kernel and ca is referred to categorical
one.

For _Health MNIST irregularly sampled_ dataset, _time_ , _id_ , and _gender_
were used as covariates for LLSM plus the intensity of the point process for
LLPPSM. When using LLSM, we employed the following structure

$f_{\text{ca}}(\text{id})+f_{\text{se}}(\text{time})+f_{\text{ca}\times\text{se}}(\text{id}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{gender}\times\text{time}),$

whereas for LLPPSM

$f_{\text{ca}}(\text{id})+f_{\text{se}}(\text{time})+f_{\text{ca}\times\text{se}}(\text{id}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{gender}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{id}\times\text{intensity})+f_{\text{ca}\times\text{se}}(\text{gender}\times\text{intensity}).$

For both variants, we used the Adam optimiser (Kingma & Ba, 2015) as
implemented in Pytorch (Paszke et al., 2019), with a learning rate equal to
0.001, which was selected based on cross-validation. After having pretrained a
standard VAE, we trained both LLSM and LLPPSM on 1000 epochs, employing early
stopping.

When training LLPPSM, we defined separate $\beta$ parameters for “healthy” and
“sick” instances and optimised them jointly together with other parameters.
For the temporal point process part of the model, the number of previous
timestamps, $D$, is 15.

###  E.2 Physionet dataset

We selected 2000 patients for training, 1917 for validation and performed
future prediction on 100 patients, not included in training and validation
sets. We used the following covariates: _time_ , _id_ , _type of ICU_ ,
_gender_ and _in-hospital mortality_ , with the corresponding additive kernel
structure for LLSM

$\displaystyle
f_{\text{ca}}(\text{id})+f_{\text{se}}(\text{time})+f_{\text{ca}\times\text{se}}(\text{id}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{gender}\times\text{time})$
$\displaystyle+f_{\text{ca}\times\text{se}}(\text{type of
ICU}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{in-hospital
mortality}\times\text{time}),$

and for LLPPSM, including intensity as a additional variable

$\displaystyle
f_{\text{ca}}(\text{id})+f_{\text{se}}(\text{time})+f_{\text{ca}\times\text{se}}(\text{id}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{gender}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{type
of ICU}\times\text{time})$ $\displaystyle+$ $\displaystyle
f_{\text{ca}\times\text{se}}(\text{in-hospital
mortality}\times\text{time})+f_{\text{ca}\times\text{se}}(\text{id}\times\text{intensity})+f_{\text{ca}\times\text{se}}(\text{gender}\times\text{intensity})$
$\displaystyle+$ $\displaystyle f_{\text{ca}\times\text{se}}(\text{type of
ICU}\times\text{intensity})+f_{\text{ca}\times\text{se}}(\text{in-hospital
mortality}\times\text{intensity}),$

The optimisation was done by Adam optimiser (Kingma & Ba, 2015) using Pytorch
(Paszke et al., 2019), with the learning rate 0.001. Both models were trained
for 400 epochs, employing early stopping, after having pretrained them by
standard VAE.

When training LLPPSM, we defined separate $\beta$ parameters based on in-
hospital mortality attribute and optimised them jointly with other parameters.
For temporal point process part, the number of previous timestamps, $D$, is
15.

Moreover, we found it useful to drop the dependence from $z^{m}$ to $y$ as was
discussed in Section 3.2. Our intuition is that in this case, if the
dependence is left, we are obliged to apply predictive distribution (Appendix
D) for both $z^{y}$ and $z^{m}$, which, due to the highly complex and sparse
nature of this dataset, cannot be modelled accurately, leading to the
accumulation of additional error.

##  Appendix F Neural network architectures

Table 9 describes neural network architecture used for both Health MNIST
setups that consists of convolutional and feedforward layers. Table 10
describes neural network architecture for the Physionet dataset. In this case,
we employed a multi layered perceptron (MLP).

Table 9: Neural Network architecture used in Health MNIST variants | Hyperparameter | Value
---|---|---
Inference network of both $z^{y}$ and $z^{m}$ | Dimensionality of input | $36\times 36$
Number of filters per convolution layer | $144$
Kernel size | $3\times 3$
Stride | $2$
Pooling | Max pooling
Pooling kernel size | $2\times 2$
Pooling stride | $2$
Number of feedforward layers | $2$
Width of feedforward layers | $300,30$
Dimensionality of latent space | $L$
Activation function of layers | RELU
Generative network of both $y$ and $m$ | Dimensionality of input | $L$
Number of transposed convolution layers | $3$
Number of filters per transposed convolution layer | $256$
Kernel size | $4\times 4$
Stride | $2$
Number of feedforward layers | $2$
Width of feedforward layers | $30,300$
Activation function of layers | RELU
Table 10: Neural Network architecture used in Physionet dataset | Hyperparameter | Value
---|---|---
Inference network of both $z^{y}$ and $z^{m}$ | Dimensionality of input | $37$
Number of feedforward layers | $2$
Width of feedforward layers | $300,30$
Dimensionality of latent space | $L$
Activation function of layers | RELU
Generative network of both $y$ and $m$ | Dimensionality of input | $L$
Number of feedforward layers | $3$
Width of feedforward layers | $30,30,300$
Activation function of layers | RELU

##  Appendix G Supplementary figures

Figure 6: Detailed overview of our model Figure 7: Inferred mean intensity
function of the point process