---
abstract: 'Electronic health records (EHR) consist of longitudinal clinical observations portrayed with sparsity, irregularity, and high-dimensionality, which become major obstacles in drawing reliable downstream clinical outcomes. Although there exist great numbers of imputation methods to tackle these issues, most of them ignore correlated features, temporal dynamics and entirely set aside the uncertainty. Since the missing value estimates involve the risk of being inaccurate, it is appropriate for the method to handle the less certain information differently than the reliable data. In that regard, we can use the uncertainties in estimating the missing values as the fidelity score to be further utilized to alleviate the risk of biased missing value estimates. In this work, we propose a novel variational-recurrent imputation network, which unifies an imputation and a prediction network by taking into account the correlated features, temporal dynamics, as well as the uncertainty. Specifically, we leverage the deep generative model in the imputation, which is based on the distribution among variables, and a recurrent imputation network to exploit the temporal relations, in conjunction with utilization of the uncertainty. We validated the effectiveness of our proposed model on two publicly available real-world EHR datasets: PhysioNet Challenge 2012 and MIMIC-III, and compared the results with other competing state-of-the-art methods in the literature.'
author:
- 'Ahmad Wisnu Mulyadi, Eunji Jun, and Heung-Il Suk,  [^1]'
bibliography:
- 'main.bib'
title: 'Uncertainty-Aware Variational-Recurrent Imputation Network for Clinical Time Series'
---

[Mulyadi : Uncertainty-Aware Variational-Recurrent Imputation Network for Clinical Time Series]{}

Bioinformatics, Deep generative model, Deep learning, Electronic health records (EHR), In-hospital mortality prediction, Missing value imputation, Time-series modeling, Uncertainty

Introduction
============

health records (EHRs) store longitudinal data consisting of patients’ clinical observations in the intensive care unit (ICU). Despite the surge of interest in clinical research on EHR, it still holds diverse challenging issues to be tackled, these include high-dimensionality, temporality, sparsity, irregularity, and bias [@cheng2016; @yadav2018]. Specifically, sequences of multidimensional medical measurements are recorded irregularly in terms of its variables and times. The reasons behind these typical measurements are diverse, such as lack of collection, documentation, or even recording faults [@wells2013; @cheng2016]. Since it carries essential information regarding a patient’s health status, improper handling of missing values might cause an unintentional bias [@wells2013; @jones2017], yielding an unreliable downstream analysis and verdict.

Complete-case analysis is an approach that draws clinical outcomes by disregarding the missing values and relying only on the observed values. However, excluding the missing data yields poor performance at high missing rates and also requires modeling separately for the distinct dataset. In fact, those missing values reflect the decisions made by health-care providers [@lipton2016]. Therefore, the missing values and their patterns contain information regarding a patient’s health status [@lipton2016] and correlate with the outcomes or target labels [@che2018]. Thus, we resort to the imputation approach by exploiting those missing patterns to improve the prediction of the clinical outcomes as the downstream task.

There exist numerous strategies for imputing missing values in the literature. Generally, imputation methods can be classified into a deterministic or stochastic approach, depending on the use of randomness [@brick1996]. Deterministic methods, such as mean [@little1987] and median filling [@acuna2004], produce only one possible value when estimating the missing values. However, it is desirable for the imputation model to generate values or samples by considering the distribution of the available observed data. Thus, it leads us to the employment of stochastic-based imputation methods.

The recent rise of deep learning models has offered potential solutions in accommodating such circumstances. Variational autoencoders (VAEs) [@kingma2014] and generative adversarial networks (GANs) [@goodfellow2014] exploit the latent distribution of high-dimensional incomplete data and generate comparable data points as the approximation estimates for the missing or corrupted values [@nazabal2018; @luo2018; @jun2019]. However, such deep generative models are insufficient for estimating the missing values of multivariate time series owing to their nature of ignoring temporal relations between a span of time points. On the other hand, by virtue of recurrent neural networks (RNNs), which have proved to perform remarkably well in modeling sequential data, we can estimate the complete data by taking into account the temporal characteristics. In this approach, [@che2018] introduced a modified gated-recurrent unit (GRU) cell to model missing patterns in the form of masking vectors and temporal delays. Likewise, BRITS [@cao2018] exploited the temporal relations of bidirectional dynamics by considering feature correlations in estimating the missing values.

Even though such models employed the stochastic approach for inferring and generating samples by utilizing both features and temporal relations, they scarcely exploited the uncertainty in estimating the missing values in multivariate time series data ([[*i.e.*]{}]{}, since the imputation estimates are not thoroughly accurate, we may introduce a fidelity score denoted by the uncertainty, which enhances the downstream task performance by emphasizing the reliable information more than the less certain information) [@he2010; @gemmeke2010; @jun2019]. We can use the imputation model, which captures the aleatoric uncertainty in estimating the missing values by placing a distribution over the output of the model [@kendall2017]. In particular, we would like to estimate the heteroscedastic aleatoric uncertainties, which are useful in cases where observation noises vary with the input [@kendall2017].

![image](img/V-RIN.pdf){width="100.00000%"}

In this work, we define our primary task as the prediction of in-hospital mortality on clinical time series data. However, since such data are portrayed with sparse and irregularly-sampled characteristics, we devise an imputation model as the secondary problem to enhance the clinical outcome predictions. We propose a novel variational-recurrent imputation network (V-RIN), as illustrated in Fig. \[fig:architecture\], which unifies the imputation and prediction networks for multivariate time series EHR data, governing both correlations among variables and temporal relations. Specifically, given the sparse data, an inference network of VAEs is employed to capture data distribution in the latent space. From this, we employ a generative network to obtain the reconstructed data as the imputation estimates for the missing values and the uncertainty indicating the imputation fidelity score. Then, we integrate the temporal and feature correlations into a combined vector and feed it into a novel *uncertainty-aware GRU* in the recurrent imputation network. Finally, we obtain the mortality prediction as a clinical verdict from the complete imputed data. In general, our main contributions in this study are as follows:

-   [We estimate the missing values by utilizing a deep generative model combined with a recurrent imputation network to capture both feature correlations and the temporal dynamics jointly, yielding the uncertainty.]{}

-   [We effectively incorporate the uncertainty with the imputation estimates in our novel uncertainty-aware GRU cell for better prediction results.]{}

-   [We evaluate the effectiveness of the proposed models by training the imputation and prediction networks jointly in an end-to-end manner, achieving superior performance on real-world multivariate time series EHR data compared to other competing state-of-the-art methods.]{}

This study extends the preliminary work published in [@jun2019]. Unlike to the preceding study, we have further expanded the proposed model by introducing more complex recurrent imputation networks, which utilize the uncertainties, instead of vanilla RNNs. We also include two additional real-world EHR datasets in our experiments, described in Section (\[section:experiments\_dataset\]), and validate the robustness of our proposed networks by comparing them with existing state-of-the-art models in the literature (Section \[section:experiments\_compmodels\]). Furthermore, we have conducted extensive experiments to discover the impacts of missing value estimation by utilizing the uncertainties when performing the downstream task (Section \[section:experiments\_ablation\] - \[section:experiments\_imputation\]).

The rest of the paper is organized as follows. In Section \[section:related\_work\], we discuss the works closely related to our proposed model in imputing the missing values. In Section \[section:methods\], we detail our proposed model. In Section \[section:experiments\], we report on the experimental results and analysis by comparing them with state-of-the-art methods. Finally, we conclude the work in Section \[section:conclusion\].

Related Work {#section:related_work}
============

Imputation strategies have been extensively devised to resolve the issue of sparse high-dimensional time series data by means of statistics, machine learning, or deep learning methods. For instance, previous works exploited statistical attributes of observed data, such as mean [@little1987] and median filling [@acuna2004], which clearly ignored the temporal relations and the correlations among variables. From the machine learning approaches, the expectation-maximization (EM) algorithm [@dempster1977], $k$-nearest neighbors (KNN) [@troyanskaya2001], and principal component analysis (PCA) [@oba2003; @mohamed2009] were proposed by considering the relationships of the features either in the original or in the latent space. Furthermore, multiple imputation by chained equations (MICE) [@white2011; @azur2011] introduced variability by means of repeating the imputation process multiple times. However, these methods ignore the temporal relations as crucial attributes in time series modeling.

The deep learning-based imputation models, are closely related to our proposed models. A previous study [@nazabal2018] leveraged VAEs to generate stochastic imputation estimates by exploiting the distribution and correlations of features in the latent space. However, it ignored the temporal relations as well as the uncertainties. Recently, GP-VAE [@fortuin2019] was proposed to obtain the latent representation by means of VAEs and model temporal dynamics in the latent space using the Gaussian process. However, since the model is merely focused on the imputation task, they required a separate model for the further downstream outcome.

To deal with time series data, a series of RNN-based imputation models were proposed. GRU-D [@che2018] considered the temporal dynamics by incorporating the missing patterns, together with the mean imputation, and forward filling with past values using temporal decay factor. Similarly, [@luo2018] trained the RNNs using such temporal decay factor and further incorporated an adversarial scheme of GANs as the stochastic approach. In the meantime, BRITS [@cao2018] were proposed to combine the feature correlations and temporal dynamic networks using bidirectional dynamics, which enhanced the accuracy by estimating missing values in both forward and backward directions. By considering the delayed gradients of the missingness in the forward and backward directions, their models were able to achieve more accurate missing values imputations. Likewise, M-RNN [@yoon2017] utilized bidirectional recurrent dynamics by operating interpolation (intra-stream) and imputation (inter-stream). Despite temporal dynamics and stochastic methods being considered in the models, the uncertainties for imputation purposes were scarcely incorporated.

As we are unsure of the actual values, we argue that these uncertainties are beneficial and can be utilized to estimate the missing values. Such uncertainty can be captured by accommodating the distribution over the model output [@kendall2017]. For this purpose, we exploited VAEs [@kingma2014] as the Bayesian networks, which are able to model the data distribution. In this work, we introduce the uncertainty as the imputation fidelity of estimates, which compensates for the potential impairment of imputation estimates. Therefore, we assumed that our model could provide reliable estimates while giving less attention to the unreliable ones determined by its uncertainties. We expect to obtain better estimates of the missing values leading to a better prediction performance as a downstream task. However, since VAEs alone are not designed to model the temporal dynamics, we combined the model with the recurrent imputation network by further utilizing those uncertainties. Thus, our proposed model differs from the aforementioned models in which that it integrates the imputation and prediction networks jointly, and the utilization of the uncertainties serve as the major motivation in our works.

Proposed Methods {#section:methods}
================

Our proposed model architecture consists of two key networks – an imputation and a prediction network – as depicted in Fig. \[fig:architecture\]. The imputation network is devised based on VAEs to capture the latent distribution of the sparse data by means of its inference network ([[*i.e.*]{}]{}, encoder $\mathcal{E}$). Then, the subsequent generative network of VAEs ([[*i.e.*]{}]{}, decoder $\mathcal{D}$) estimates the reconstructed data distribution. We regard its reconstructed values as the imputation estimates while exploiting its variances as the uncertainties to be further utilized in the recurrent imputation network sequentially.

The succeeding recurrent imputation networks are built upon RNNs to model the temporal dynamics. For each time step, we employ the regression layer to model the temporal relations incorporated within the hidden states of RNNs cell in imputing the missing values. However, as we also consider the time gap between observed values, we incorporate the time decay factor, which is then exploited in such hidden states, leading to *decayed hidden states*. Eventually, by systematically unifying the imputation estimates obtained from VAEs and RNNs, we expect to acquire more likely estimates by considering feature correlations and temporal relations over time, including the utilization of the uncertainty. By doing so, we expect to acquire more reliable prediction outcomes. We describe each of the networks more specifically in the following sections after introducing the data representation.

Data Representation
-------------------

Given the multivariate time series EHR data of $N$ number of patients, a set of clinical observations and their corresponding label is denoted as $\{\rmX^{(n)}, y^{(n)}\}^{N}_{n=1}$. For each patient, we have $\rmX^{(n)} = [\rvx^{(n)}_1, \ldots, \rvx^{(n)}_t, \ldots, \rvx^{(n)}_T]^\top \withdim{T \times D}$, where $T$ and $D$ represent the time points and variables, respectively; $\rvx^{(n)}_t \withdim{D}$ denotes all observable variables at time point $t$, and $x^{d,(n)}_t$ is the $d$-th element of variables at time point $t$. In addition, it has the corresponding clinical label $y^{(n)} \in \{0,1\}$, representing the clinical outcomes. In our case, it denotes the in-hospital mortality of a patient, which falls into a binary classification problem. For the sake of clarity, we omit the superscript $(n)$ hereafter.

As $\rmX$ is characterized with sparsity properties, we address the missing values by introducing masking matrix $\rmM \in \{0,1\}^{T\times D}$, indicating whether values are observed or missing. In addition, we define a new data representation $\tilde{\rmX} = [\tilde{\rvx}_1, \ldots, \tilde{\rvx}_t, \ldots, \tilde{\rvx}_T]^\top \withdim{T \times D}$ to be fed into the model, where we initialize the missing value with zero [@lipton2016; @nazabal2018; @jun2019] as follows: $$\begin{aligned}
    m^d_t &= 
    \begin{cases}
        1 &\quad \textrm{if} \ x^d_t \ \textrm{is observed},\\
        0 &\quad \textrm{otherwise}
    \end{cases},
    \\
    \tilde{x}^d_t &= 
    \begin{cases}
        x^d_t &\quad \textrm{if} \ m^d_t = 1,\\
        0 &\quad \textrm{otherwise.}
    \end{cases}\end{aligned}$$ Another consideration in dealing with the sparse data is that there exists a time gap between two observed values. Such information in fact carries a piece of essential information in estimating the missing values over the times. Thus, we accommodate this information by further devising the time delay matrix ${\mDelta} \withdim{T\times D}$, which is derived from $\rvs \withdim{T}$, denoting the timestamp of the measurement. As initialization, we fix ${\vdelta}_t = \mOne$ for the $t=1$, while setting the time delay for the remaining times $(t > 1)$ by referring to a masking matrix and a timestamp vector as follows: $$\begin{aligned}
    \Delta^d_t = 
    \begin{cases}
        s_t - s_{t-1} &\quad \textrm{if} \ m^d_{t-1}=1, \\
        s_t - s_{t-1} + \Delta^d_{t-1} &\quad \textrm{otherwise.} \\
    \end{cases}\end{aligned}$$

![image](img/GRUs.pdf){width="0.8\linewidth"}

VAE-based Imputation Network
----------------------------

Given the observations at each time point $\tilde{\rvx}_t$, we infer $\rvz\withdim{k}$ as the latent representation with $k$ as its corresponding dimension by making use of the inference network, using the true posterior distribution $p_\phi(\rvz|\tilde{\rvx}_t)$. Intuitively, we assume that $\tilde{\rvx}_t$ is generated from some unobserved random variable $\rvz$ by some conditional distribution $p_\theta(\tilde{\rvx}_t|\rvz)$, while $\rvz$ is generated from a prior distribution $p_\theta(\rvz)$, which can be interpreted as the hidden health status of the patient. In addition, we define the marginal likelihood as $p_\theta(\tilde{\rvx}_t)$ by integrating out the joint distribution $p_\theta(\tilde{\rvx}_t,\rvz)$ for $\rvz$ defined as

$$p_\theta(\tilde{\rvx}_t) = \int p_\theta(\tilde{\rvx}_t,\rvz)d\rvz = \int p_\theta(\tilde{\rvx}_t|\rvz) p_\theta(\rvz) d\rvz.$$

However, in practice, this is analytically intractable. Since, $p_\theta(\rvz|\tilde{\rvx}_t) = p_\theta(\tilde{\rvx}_t|\rvz)p_\theta(\rvz)/p_\theta(\tilde{\rvx}_t)$, the true posterior becomes intractable as well. Therefore, we approximate it with $q_\phi(\rvz|\tilde{\rvx}_t)$ using a Gaussian distribution $\mathcal{N}(\vmu_z, \textrm{diag}(\vsigma^2_z))$, where the mean and log-variance are obtained such that $$\begin{aligned}
    \vmu_z &= \gE_\mu(\tilde{\rvx}_t;\phi), \quad \log \vsigma^2_z = \gE_\sigma(\tilde{\rvx}_t;\phi),\end{aligned}$$ where $\gE_{\{\mu, \sigma\}}$ denotes the inference network with parameter $\phi$. Furthermore, we apply the reparameterization trick [@kingma2014] as $\rvz = \vmu_z + \vsigma_z \odot \ \vepsilon_z$, where $\vepsilon_z \sim \mathcal{N}(\mZero,\mIdentity)$, with $\odot$ denoting element-wise multiplication, thus, making it possible to be differentiated and trained using standard gradient methods.

Furthermore, given this latent vector $\rvz$, we estimate $p_\theta(\tilde{\rvx}_t|\rvz)$ by means of the generative network $\gD$ with parameters $\theta$ as $$\begin{aligned}
    \vmu_{x,t} &= \gD_\mu(\rvz;\theta), \quad \log \vsigma^2_{x,t} = \gD_\sigma(\rvz;\theta),\end{aligned}$$ where $\vmu_{x,t}$ and $\vsigma^2_{x,t}$ denote the mean and variance of reconstructed data distribution, respectively. We regard the mean as the estimate to the missing values and maintain the observed values in $\bar{\rvx}_t$ by making use of corresponding masks vector as $$\begin{aligned}
    \bar{\rvx}_t = \rvm_t \odot \tilde{\rvx}_t + (\mOne - \rvm_t) \odot \vmu_{x,t}.\end{aligned}$$ In the meantime, we regard the variance of reconstructed data as the uncertainty to be further utilized in the recurrent imputation process. For this purpose, we introduce an uncertainty matrix $\bar{\rmU} \withdim{T\times D}$ with $\hat{\mSigma}_x= [\textrm{diag}(\vsigma_{x,1}), \ldots, \textrm{diag}(\vsigma_{x,t}), \ldots, \textrm{diag}(\vsigma_{x,T})]^\top \withdim{T \times D}$. We quantify this uncertainty as the fidelity score of the missing value estimates. In particular, we set the corresponding uncertainty to zero if the data is observed, indicating that we are confident with full trust in the observation, and set it as a value $\sigma^d_{x,t}$ if the corresponding value is missing as $$\begin{aligned}
    \bar{\rvu}_t = (\mOne - \rvm_t) \odot \hat{\vsigma}_{x,t}\end{aligned}$$

Finally, as an output of this VAE-based imputation network, we acquire the set $\{\bar{\rmX}, \bar{\rmU}\}$ denoting the imputed values and their corresponding uncertainty, respectively. Furthermore, to alleviate the bias of missing value estimations, we utilize this uncertainty in the following recurrent imputation network.

Recurrent Imputation Network
----------------------------

The recurrent imputation network is based on RNNs, where we further model the temporal relations in the imputed data and exploit the uncertainties. While both GRU [@cho2014] (depicted in Fig. \[fig:grus\_cells\]a) and long-short term memory (LSTM) [@hochreiter1997] are feasible choices, inspired by the previous work of [@che2018] depicted in Fig. \[fig:grus\_cells\]b, we employed a modified GRU cell by leveraging the uncertainty-aware GRU cell to consider further uncertainty and the temporal decaying factor, which is depicted in Fig. \[fig:grus\_cells\]c. Graphical illustrations of several forms of GRU cells are depicted in Fig. \[fig:grus\_cells\].

Specifically, at each time step $t$, we produce the uncertainty decay factor $\vupsilon_t$ in Eq. (\[eq:unc\_decay\]) using a negative exponential rectifier to guarantee $\vupsilon_t \in (\mZero,\mOne]$ [@che2018; @cao2018]. $$\label{eq:unc_decay}
    \vupsilon_t = \textrm{exp}\{-\textrm{max}(\mZero, \rmW_u \bar{\rvu}_t + \rvb_u)\}$$

We utilize such factors to emphasize the reliable estimates and give less attention to the uncertain ones. In particular, we first employ a fully-connected layer to $\bar{\rvx}_t$ and element-wise multiply with the uncertainty decay factor $\vupsilon_t$ as follows:

$$\label{eq:unc_decay_odot_x}
    \rvx^\upsilon_t = (\rmW_{\upsilon}\bar{\rvx}_t + \rvb_{\upsilon}) \odot \vupsilon_t.$$

Note that we zeroed-out the diagonal of the parameter $\rmW_{\upsilon}$ to enforce the estimations based on other features. Thus, we obtain $\rvx^\upsilon_t$ as the *feature-based correlated estimates* to the missing values.

In addition, we further consider the missing value estimates based on the temporal relations. For this purpose, we employ the time delay $\vdelta_t$ which is an essential element to capture temporal relations and missing patterns of the data [@che2018]. We exploit this information as the temporal decay factor $\vgamma_t \in (\mZero,\mOne]$ as follows $$\label{eq:time_decay}
    \vgamma_t = \textrm{exp}\{-\textrm{max}(\mZero, \rmW_{\gamma} \vdelta_t + \rvb_{\gamma})\}.$$ Meanwhile, by employing the GRU, we obtain the hidden state $\rvh$ as the comprehensive information compiled from the preceding sequences. Thus, we take advantage of the temporal decay factor in governing the influence of past observations embedded into hidden states using the form of decayed hidden states as $$\label{eq:decayed_hidden_states}
    \hat{\rvh}_{t-1} = \rvh_{t-1} \odot \vgamma_t.$$ Thereby, given the previous hidden states $\hat{\rvh}_{t-1}$, we can estimate the current complete observation $\rvx^r_t$ through regression. $$\label{eq:x_reg}
    \rvx^r_t = \rmW_r\hat{\rvh}_{t-1} + \rvb_r.$$ In addition, we further make use of those estimates by applying another operation:

$$\label{eq:x_temporal}
    \rvx^{\tau}_t = \rmW_{\tau}\rvx^r_t + \rvb_{\tau},$$

again setting the diagonal parameter of $\rmW_{\tau}$ to be zeros to consider *the feature-based estimation on top of the temporal relations from the previous hidden states*.

Hence, we have a pair of imputed values $\{\rvx^\upsilon_t, \rvx^{\tau}_t \}$, corresponding to missing value estimates obtained from the VAE by considering the uncertainties, and from the recurrent imputation network, respectively. We then merge this information jointly to get combined vector $\rvc_t$ comprising both estimates by simply employing a $1\times1$ convolution operation ($\ast$) as $$\label{eq:x_comb}
    \rvc_t = \rvx^{\upsilon}_t \ast \rvx^{\tau}_t$$ Finally, we obtain the complete vector $\rvx^c_t$ by replacing the missing values with the combined estimates as $$\label{eq:x_final}
    \rvx^c_t = \rvm_t \odot \tilde{\rvx}_t + (\mOne - \rvm_t) \odot \rvc_t.$$ In addition, we concatenate the complete vector with the corresponding mask, and then feed it into the modified GRU cell to obtain the subsequent hidden states $$\label{eq:rnn}
    \rvh_t = \sigma(\rmW_h\hat{\rvh}_{t-1} + \rmV_h[\rvx^c_t \circ \rvm_t] + \rvb_h).$$ Lastly, to predict the in-hospital mortality as the clinical outcome, we utilize the last hidden state $\rvh_T$ to get the predicted label $\hat{y}$ such that $$\label{eq:output}
    \hat{y} = \sigma(\rmW_y\rvh_T + \rvb_y).$$ Hereby, $\rmW_{\{u, \upsilon, \gamma, r, \tau, h, y\}}$, $\rmV_h$, and $\rvb_{\{u, \upsilon, \gamma, r, \tau, h, y\}}$ are the learnable parameters in our recurrent imputation network.

$\nabla_{\vTheta} \leftarrow \mZero$\

Learning
--------

We describe the composite loss function, comprising the imputation and prediction loss function to tune all model parameters jointly, which are $\vTheta = \{\theta, \phi, \rmW_{\{u, \upsilon, \gamma, r, \tau, h, y\}}, \rmV_h, \rvb_{\{u, \upsilon, \gamma, r, \tau, h, y\}}\}$. Such loss function accommodates the VAEs and the recurrent imputation network as well. By means of VAEs, we define the loss function $\gL_{vae}$ to maximize the variational evidence lower bound (ELBO) that comprises the reconstruction loss term and the Kullback-Leibler divergence. We add $\ell_1$-regularization to introduce sparsity into the network with $\lambda_1$ as the hyperparameter. Moreover, for each time step, we measure the difference between the observed data and the combined imputation estimates by the mean absolute error (MAE) as $\gL_{reg}$. $$\label{eq:loss_vae}
    \gL_{vae} =  \sum^N_{n=1} \sum^T_{t=1} \mathbb{E}_{q_\phi(\rvz|\tilde{\rvx}^{(n)}_t)}[\log p_\theta(\tilde{\rvx}^{(n)}_t|\rvz)]$$ $$-  \KLDiv{q_\phi(\rvz|\tilde{\rvx}^{(n)}_t)}{p_\theta(\rvz)} + \lambda_1 \lVert(\theta;\phi)\lVert_1$$ $$\label{eq:loss_reg}
    \gL_{reg} = \sum^N_{n=1}{\sum^T_{t=1}{ \gL_{MAE}(\tilde{\rvx}^{(n)}_t \odot \rvm^{(n)}_t, \rvc^{(n)}_t \odot \rvm^{(n)}_t)}}$$ Furthermore, we define the binary cross-entropy loss function $\gL_{pred}$ to evaluate the prediction of in-hospital mortality. Thus, we define the overall composite loss function $\gL_{total}$ as $$\label{eq:loss_total}
    \gL_{total} = \alpha \ \gL_{vae} + \beta \ \gL_{reg} + \gL_{pred},$$ where $\alpha$ and $\beta$ are the hyperparameters to represent the ratio between the $\gL_{vae}$ and $\gL_{reg}$, respectively.

Note that our proposed model is also applicable to consider the bidirectional dynamics. Such a scenario can be carried out by having the forward and backward direction of the data fed into the recurrent imputation network. By doing so, we make our proposed model a fair comparison to M-RNN [@yoon2017], BRITS-I, and BRITS [@cao2018] which adopted such strategy to achieve better estimates of the missing values and the prediction outcomes. In bidirectional cases, we employ an additional consistency loss function to the aforementioned total loss function $\gL_{total}$, as introduced in [@cao2018] to impose consistency estimates for each time step in both directions as $$\label{eq:loss_consistency}
    \gL_{cons} = \sum^N_{n=1}{\sum^T_{t=1}{ \gL_{MAE}(\rvx^{c,(n)}_t, \rvx^{c',(n)}_t)}} ,$$ with $\rvx^{c,(n)}_t$ and $\rvx^{c',(n)}_t$ denoting the estimates from the forward and backward direction, respectively. Another hyperparameter $\xi$ could be introduced for this consistency loss to optimize the model.

Lastly, we use stochastic gradient descent in an end-to-end manner to optimize the model parameters during the training. We summarize the overall training steps of our proposed framework in Algorithm \[algo:model\].

Experiments {#section:experiments}
===========

Dataset and Implementation Setup {#section:experiments_dataset}
--------------------------------

### **PhysioNet 2012 Challenge**

PhysioNet[^2] [@goldberger2000; @ikaro2012] consists of 35 irregularly sampled clinical variables ([[*e.g.*]{}]{},  heart and respiration rate, blood pressure, [[*etc.*]{}]{}) from 4,000 patients during their first 48 hours of medical care in the ICU. Note that we ignore the demographic information and categorical data types from this dataset. Hereby, we exploit only the clinical time series data. From those samples, we excluded three patients with no observations at all. We sampled the observations hourly, using the time window as the timestamps, and took the average of values in cases of multiple measurements within this time window. It resulted in sparse EHR data with an average missing rate of 80.51%. Our aim is to predict the in-hospital mortality of patients, with 554 positive mortality labels (13.86%). As for the implementation setup of the PhysioNet dataset, we employed three layers of feedforward networks for the inference network of VAEs with hidden units of $\{64,24,10\}$, with $10$ denoting the dimension of latent representation. The generative network has equal numbers of hidden units with those of the inference network but in the reverse order. We employed hyperbolic tangent (tanh) as the non-linear activation function for each hidden layer. Prior to those activation functions, we also applied batch normalization and a dropout rate of $\{0.1,0.3\}$ for classification and imputation tasks, respectively. We employed modified GRU for recurrent imputation network with $64$ hidden units.

### **MIMIC-III**

The Medical Information Mart for Intensive Care (MIMIC-III)[^3] [@johnson2016] dataset consists of 53,432 ICU stays of adult patients in the Beth Israel Deaconess Medical Center in the period of 2001 – 2012 [@johnson2016]. We selected 99 variables from several source tables, such as laboratory tests, inputs to patients, outputs from patients, and drug prescriptions tables, resulting in a cohort of 13,998 patients with 1,181 positive in-hospital mortality labels (8.44%) among them. Moreover, from those irregular measurements, we further sampled the data into two hourly samples for the first 48 hours of their medical care, leading to an average missing rate of 93.92%. As in the case of PhysioNet, we took the average value if there existed multiple measurements. We referred to [@che2018; @purushotham2018] in pre-processing the MIMIC-III cohort. We employed three layers of feedforward networks for the inference network of VAEs with hidden units of $\{128,32,16\}$ and an equal number of hidden units in reverse for the generative network. Likewise, for each hidden layer, we used batch normalization, drop out rate of $0.3$ for both of classification and imputation tasks, and tanh activation function. Furthermore, $64$ hidden units were employed for the recurrent imputation network.

We trained the proposed model on both datasets using an Adam optimizer with $100$ epochs and $64$ mini-batches. For imputation task, we fixed the learning rates of $0.0005$ and $0.0003$ for PhysioNet and MIMIC-III, respectively, while $0.005$ for both datasets on the classification task. We set $\lambda_{1}$ and weight decay equally with $1e-5$. For bidirectional models, we set $\xi=0.1$ as the hyperparameter of the consistency loss.

Tasks {#section:experiments_task}
-----

In this work, we validated the performance of our proposed models from two perspectives: (1) the in-hospital mortality prediction (classification), and (2) missing value imputation.

### Classification Task

Our primary goal for this work is to predict the in-hospital mortality as the binary classification task. For this purpose, we reported the test result on the in-hospital mortality prediction task from the 5-fold cross-validation in terms of the average of the Area Under the Curve (AUC) ROC. Additionally, to measure the robustness of the models in dealing with the imbalance data portrayed in both datasets, we also reported the Area Under the Precision-Recall Curve (AUPRC). We randomly removed samples of the training data with $\{5\%, 10\%\}$ scenarios and left the validation and test set untouched to be reported as the results.

### Imputation Task

For the secondary task, we additionally evaluated imputation performance on the missing values. As for this task, we randomly removed samples with settings of $\{5\%, 10\%\}$ of observed data from all training, validation, and test set as the ground truth. Then we reported the test result from the 5-fold cross-validation by measuring the MAE. In addition, we also measured the other most-used imputation similarity metric, which is mean relative error (MRE). Given $\rvx^{(i)}$ and $\rvx^{c,(i)}$ as the ground truth and imputation estimates of $i$-th item, respectively, and $N_{GT}$ ground truth items in total, MAE and MRE defined as $$MAE = \frac{\sum_{i=1}^{N_{GT}} |\rvx^{c,(i)} - \rvx^{(i)}|}{N_{GT}},$$ $$MRE = \frac{\sum_{i=1}^{N_{GT}} |\rvx^{c,(i)} - \rvx^{(i)}|}{\sum_{i=1}^{N_{GT}} |\rvx^{(i)}|}.$$

![image](img/ablation_studies_physionet.pdf){width="95.00000%"}

  -- -- -- -- -- -- --
                    
                    
                    
                    
                    
                    
                    
                    
  -- -- -- -- -- -- --

Comparative Models {#section:experiments_compmodels}
------------------

We compared the performance of our proposed model in carrying the aforementioned tasks with the closely-related competing state-of-the-art models in the literature by grouping them into unidirectional and bidirectional models.

### Unidirectional Models

-   [[**GRU-D**]{} [@che2018] estimates the missing values by utilizing the informative missing value patterns in the form of the masking and time decay factor using modified GRU cells.]{}

-   [[**RITS-I**]{} [@cao2018] utilizes unidirectional dynamics that rely solely on temporal relations through regression operations.]{}

-   [[**RITS**]{} [@cao2018] is devised based on RITS-I by further taking into account the feature correlations function. Furthermore, it utilizes the temporal decay as the factor to weigh between both features and temporal-based estimates.]{}

-   [[**V-RIN**]{} (Ours) is based on our proposed model except that we ignored the uncertainty decay factor. Specifically, we excluded Eq. (\[eq:unc\_decay\]) and omitted the element-wise multiplication operation of the $\vupsilon_t$ in Eq. (\[eq:unc\_decay\_odot\_x\]).]{}

-   [[**V-RIN-full**]{} (Ours) executed all operations in the proposed model including feature-based correlations, temporal relations, and the uncertainty decay utilization. ]{}

### Bidirectional Models

-   [[**M-RNN**]{} [@yoon2017] exploits the multi-directional RNNs which execute both interpretation and imputation to infer the missing data.]{}

-   [[**BRITS-I**]{} [@cao2018] is based on the RITS-I by extending it to be able to handle bidirectional dynamics in estimating the missing values.]{}

-   [[**BRITS**]{} [@cao2018] takes the bidirectional dynamics of RITS in handling the sparsity in the data. Both and BRITS additionally employ consistency loss of forward and backward directions in their attempt to estimate the missing values more precisely.]{}

-   [To make a fair comparison, we extended our proposed model of [**V-RIN**]{} and [**V-RIN-full**]{} into the bidirectional models by means of the recurrent imputation networks. Similar to BRITS-I and BRITS, we further computed consistency loss in Eq. (\[eq:loss\_consistency\]) for both forward and backward directions of the estimates. ]{}

  -- -- --------- ----------- --------- -----------
                                        
         **AUC**   **AUPRC**   **AUC**   **AUPRC**
                                        
                                        
                                        
                                        
                                        
                                        
                                        
                                        
                                        
                                        
  -- -- --------- ----------- --------- -----------

  -- -- --------- ----------- --------- -----------
                                        
         **AUC**   **AUPRC**   **AUC**   **AUPRC**
                                        
                                        
                                        
                                        
                                        
                                        
                                        
                                        
                                        
                                        
  -- -- --------- ----------- --------- -----------

Experimental Results : Ablation Studies {#section:experiments_ablation}
---------------------------------------

As part of the ablation studies, we report the performance of the unidirectional model of the V-RIN and V-RIN-full models on the in-hospital mortality classification task. Firstly, we investigated the effect of the pair $\{\alpha, \beta\}$ as the hyperparameter in Eq. (\[eq:loss\_total\]). We reflected these parameters as a ratio to weigh the imputation by the VAEs and the recurrent imputation network to achieve the optimal performance on estimating the missing values and classifying the clinical outcomes. For each parameter, we defined set of values in the range of $[0.1, 1.0]$.

For PhysioNet, as illustrated in Fig. \[fig:ablation\_study\_physionet\], in general, we observed that in almost all the combination settings, V-RIN-full achieved higher performance than V-RIN. We interpreted these findings that introducing the uncertainty helps the model in estimating the missing values, leading to better classifying the outcome. Both models were able to achieve high performance in terms of their average AUC scores of $0.8311 \pm 0.0079$ for V-RIN and $0.8347 \pm 0.0125$ for V-RIN-full. These AUC results were obtained with settings of $\{\alpha=0.75,\ \beta=0.25\}$ and $\{\alpha=1.0,\ \beta=0.75\}$, respectively. For these results, we observed that the model favored the emphasis on the feature correlations over the temporal relations to obtain its best performance. For the case of V-RIN, once we tried to increase the $\beta>0.25$, we observed that the classification performance was degraded to some degree. In contrast, the performance of V-RIN-full was considerably better when we increased the $\beta$ parameter. We also carried out similar ablation studies on the MIMIC-III dataset, which is reported in the supplementary material. To summarize, for both datasets, we argue that both features and temporal relations are essential in estimating the missing values with some latent proportion.

Table \[tab:ablation\_studies\] presents the comparison of our model with closely related models on both classification and imputation tasks, such as VRNN [@chung2015], which integrates VAEs for each time step of RNNs, and VAE+RNN [@jun2019], which employs VAEs followed by RNNs without incorporating the uncertainty. For the case of VRNN on PhysioNet and MIMIC-III, we observed that the classification performance was the lowest among reported models in terms of both its AUC and AUPRC. In addition, their imputation performance was high in comparison to both our proposed models. As in the case of VAE+RNN, in comparison to VRNN, we noticed a considerable improvement in the performance on PhysioNet and MIMIC-III, especially in terms of AUPRC. VAE+RNN is, in fact, the model closely related to ours that it executes the imputation process by first exploiting the feature correlations followed by temporal dynamics in exact order. However, [@jun2019] employed the vanilla RNNs instead of recurrent imputation network, which is a novel extension in this study. By introducing the temporal decay in V-RIN, the model was better able to learn the temporal dynamics effectively, resulting in better AUC and AUPRC, as well as better imputation results in terms of MAE and MRE on both datasets, by a large margin. Finally, once we introduced the uncertainty which is incorporated in the recurrent imputation network of V-RIN-full, we observed a significant enhancement of AUC, AUPRC, MAE, and MRE on both datasets. Thus, we conclude that the utilization of both temporal decay and the uncertainties are beneficial in both imputation and classification tasks.

  -- -- --------- --------- --------- ---------
                                      
         **MAE**   **MRE**   **MAE**   **MRE**
                                      
                                      
                                      
                                      
                                      
                                      
                                      
                                      
                                      
                                      
  -- -- --------- --------- --------- ---------

  -- -- --------- --------- --------- ---------
                                      
         **MAE**   **MRE**   **MAE**   **MRE**
                                      
                                      
                                      
                                      
                                      
                                      
                                      
                                      
                                      
                                      
  -- -- --------- --------- --------- ---------

Classification Result Analysis {#section:experiments_classification}
------------------------------

We presented the experimental results of the in-hospital mortality prediction in comparison with other competing methods in terms of average AUC and AUPRC in Table \[tab:classifcation\_physionet\] and Table \[tab:classification\_mimic\] for PhysioNet and MIMIC-III, respectively. We reported the evaluation of both unidirectional and bidirectional models for 10% and 5% removal on both datasets. For the case of unidirectional models on PhysioNet with 10% and 5% removal, our V-RIN model achieved better performance in terms of AUC and AUPRC by a large margin compared to all comparative models even with the RITS, which utilized both the features and temporal relations, sequentially. The highest AUC was obtained by the V-RIN-full in both removal percentage scenarios. As for bidirectional models, all models in the 10% and 5% scenarios improved their performance, except for the M-RNN, which achieved lowest results among all competing methods. Although M-RNN employed similar strategies using bidirectional dynamics, it struggles to perform the task properly. Both AUC and AUPRC of M-RNN was lower than GRU-D, despite GRU-D using only the forward directional data. The highest AUC was achieved by our V-RIN-full with $0.8401 \pm 0.0100$ and $0.8422 \pm 0.0070$ for the 10% and 5% removal scenarios, respectively. These findings reassure that the utilization of the uncertainty is truly beneficial in estimating the missing values.

As for MIMIC-III, we reported the classification performance in Table \[tab:classification\_mimic\]. In this dataset, we observed quite similar patterns to PhysioNet, with our model outperforming all competing models. The imbalance ratio, missing ratios, and dimensionality of MIMIC-III are much greater than PhysioNet. We found that the performance of V-RIN and V-RIN-full was comparable to some extent. In unidirectional models with 10%, our V-RIN-full achieved the highest AUC and AUPRC, whereas the V-RIN achieved the highest performance by a small margin in the bidirectional models. This contrasts with the 5% scenarios, where V-RIN-full with bidirectional models achieved the highest AUC and AUPRC of $0.8705 \pm 0.0148$ and $0.4199 \pm 0.016$, respectively.

![image](img/lrp_timeseries.pdf){width="100.00000%"}

Imputation Results Analysis {#section:experiments_imputation}
---------------------------

As the secondary task, we further evaluated the imputation performance of our proposed model in contrast with the comparative models. Table \[tab:imputation\_physionet\] presents the experimental results on the PhysioNet dataset. GRU-D which exploited the mean and temporal relations in imputing the values struggle to achieve good performance. By using bidirectional dynamics, M-RNN obtained better results than GRU-D. However, those models are still inferior to the RITS variants in both directional scenarios. Overall, both of our proposed models revealed its imputation robustness on this dataset, by exhibiting best performance, in both unidirectional and bidirectional scenario, consistently. Table \[tab:imputation\_mimic\] presents the imputation performance on MIMIC-III. In general, the results showed that our proposed model performed better than comparative models, and comparable to BRITS in bidirectional models for both 10% and 5% removal scenarios.

We further exploited layer-wise relevance propagation (LRP) [@montavon2018] through a publicly available toolbox: [^4] [@alber2018] to examine the relevance of the features. Specifically, we employed the LRP-$\alpha_1\beta_0$ rule [@montavon2018; @alber2018] to discover which features induce the activation of the neurons in carrying out the classification task. As a result, we depicted the imputed values as the input to the recurrent imputation network and further revealed the positive and negative relevance of each feature over time on both datasets in Fig. \[fig:lrp\_timeseries\]. The observed values are illustrated as black circles. By means of VAEs, we obtained the imputation estimates, which are highlighted by the hollow diamond markers alongside the uncertainties, illustrated with the shaded areas. Then, by employing the recurrent imputation network, we acquired the final estimates to the missing values depicted as the hollow circles. Overall, we observed that the relevances were getting stronger toward the end of the time period, regardless of its sign. Furthermore, compared to the observed values, the missing value estimates with high uncertainties tend to hold a low relevances. Thus, this demonstrates the benefits of the uncertainty utilization in the recurrent imputation networks for the downstream task.

Conclusion {#section:conclusion}
==========

In this study, we proposed a novel unified framework consisting of imputation and prediction networks for sparse high-dimensional multivariate time series. It combined a deep generative model with a recurrent model to capture feature correlations and temporal relations for estimating the missing values by taking into account the uncertainty. We utilized the uncertainties as the fidelity of our estimation and incorporated them for predicting clinical outcomes. We evaluated the effectiveness of the proposed model with the PhysioNet 2012 Challenge and MIMIC-III datasets as the real-world EHR multivariate time series data, proving the superiority of our model in the in-mortality prediction task, compared to other comparative state-of-the-art models in the literature.

Acknowledgment {#acknowledgment .unnumbered}
==============

This work was supported by Institute of Information & communications Technology Planning & Evaluation (IITP) grant funded by the Korea government (MSIT) (No. 2017-0- 01779, A machine learning and statistical inference framework for explainable artificial intelligence, and No. 2019-0-00079, of Artificial Intelligence (Korea University)).

[^1]: A. W. Mulyadi and E. Jun are with the Department of Brain and Cognitive Engineering, Korea University, Seoul 02841, Korea (e-mail: wisnumulyadi@korea.ac.kr, ejjun92@korea.ac.kr). H.-I. Suk is with the Department of Artificial Intelligence and the Department of Brain and Cognitive Engineering, Korea University, Seoul 02841 Korea (e-mail: hisuk@korea.ac.kr). (*Corresponding author: Heung-Il Suk*)

[^2]: Publicly available on <https://physionet.org/content/challenge-2012/1.0.0/>

[^3]: Publicly available on <https://physionet.org/content/mimiciii/1.4/>

[^4]: <https://github.com/albermax/innvestigate/>
