Title: RECURRENT NEURAL NETWORKS FOR MULTIVARI-ATE TIME SERIES WITH MISSING VALUES

Abstract: Multivariate time series data in practical applications, such as health care, geoscience, and biology, are characterized by a variety of missing values. In time series prediction and other related tasks, it has been noted that missing values and their missing patterns are often correlated with the target labels, a.k.a., informative missingness. There is very limited work on exploiting the missing patterns for effective imputation and improving prediction performance. In this paper, we develop novel deep learning models, namely GRU-D, as one of the early attempts. GRU-D is based on Gated Recurrent Unit (GRU), a state-of-the-art recurrent neural network. It takes two representations of missing patterns, i.e., masking and time interval, and effectively incorporates them into a deep model architecture so that it not only captures the long-term temporal dependencies in time series, but also utilizes the missing patterns to achieve better prediction results. Experiments of time series classification tasks on real-world clinical datasets (MIMIC-III, PhysioNet) and synthetic datasets demonstrate that our models achieve state-of-the-art performance and provides useful insights for better understanding and utilization of missing values in time series analysis.

1. INTRODUCTION:

Multivariate time series data are ubiquitous in many practical applications ranging from health care, geoscience, astronomy, to biology and others. They often inevitably carry missing observations due to various reasons, such as medical events, saving costs, anomalies, inconvenience and so on. It has been noted that these missing values are usually informative missingness (Rubin, 1976) , i.e., the missing values and patterns provide rich information about target labels in supervised learning tasks (e.g, time series classification). To illustrate this idea, we show some examples from MIMIC-III, a real world health care dataset in Figure 1 . We plot the Pearson correlation coefficient between variable missing rates, which indicates how often the variable is missing in the time series, and the labels of our interests such as mortality and ICD-9 diagnoses. We observe that the missing rate is correlated with the labels, and the missing rates with low rate values are usually highly (either positive or negative) correlated with the labels. These findings demonstrate the usefulness of missingness patterns in solving a prediction task.

In the past decades, various approaches have been developed to address missing values in time series (Schafer & Graham, 2002) . A simple solution is to omit the missing data and to perform analysis only on the observed data. A variety of methods have been developed to fill in the missing values, such as smoothing or interpolation (Kreindler & Lumsden, 2012) , spectral analysis (Mondal & Percival, 2010) , kernel methods (Rehfeld et al., 2011) , multiple imputation (White et al., 2011) , and EM algorithm (García-Laencina et al., 2010) . Schafer & Graham (2002) and references therein provide excellent reviews on related solutions. However, these solutions often result in a twostep process where imputations are disparate from prediction models and missing patterns are not effectively explored, thus leading to suboptimal analyses and predictions (Wells et al., 2013) .

In the meantime, Recurrent Neural Networks (RNNs), such as Long Short-Term Memory (LST-M) (Hochreiter & Schmidhuber, 1997) and Gated Recurrent Unit (GRU) (Cho et al., 2014) , have shown to achieve the state-of-the-art results in many applications with time series or sequential data, including machine translation (Bahdanau et al., 2014; Sutskever et al., 2014) and speech recognition (Hinton et al., 2012) . RNNs enjoy several nice properties such as strong prediction performance as well as the ability to capture long-term temporal dependencies and variable-length observations. RNNs for missing data has been studied in earlier works (Bengio & Gingras, 1996; Tresp & Briegel, 1998; Parveen & Green, 2001) and applied for speech recognition and blood-glucose prediction.

Recent works (Lipton et al., 2016; Choi et al., 2015) tried to handle missingness in RNNs by concatenating missing entries or timestamps with the input or performing simple imputations. However, there have not been works which model missing patterns into a systematically modified RNN structure for time series classification problems. Exploiting the power of customized RNN models along with the informativeness of missing patterns is a new promising venue to effectively model multivariate time series and is the main motivation behind our work.

In this paper, we develop a novel deep learning model based on GRU, namely GRU-D, to effectively exploit two representations of informative missingness patterns, i.e., masking and time interval.

Masking informs the model which inputs are observed (or missing), while time interval encapsulates the input observation patterns. Our model captures the observations and their dependencies by applying masking and time interval (using a decay term) to the inputs and network states of GRU, and jointly train all model components using back-propagation. Thus, our model not only captures the long-term temporal dependencies of time series observations but also utilizes the missing patterns to improve the prediction results. Empirical experiments on real-world clinical datasets as well as synthetic datasets demonstrate that our proposed model outperforms strong deep learning models built on GRU with imputation as well as other strong baselines. These experiments show that our proposed method is suitable for many time series classification problems with missing data, and in particular is readily applicable to the predictive tasks in emerging health care applications. Moreover, our method provides useful insights into more general research challenges of time series analysis with missing data beyond classification tasks, including 1) a general deep learning framework to handle time series with missing data, 2) effective solutions to characterize the missing patterns of not missing-completely-at-random time series data such as modeling masking and time interval, and 3) an insightful approach to study the impact of variable missingness on the prediction labels by decay analysis.

2. RNN MODELS FOR TIME SERIES WITH MISSING VARIABLES:

We denote a multivariate time series with D variables of length T as X = (x 1 , x 2 , . . . , x T ) T ∈ R T ×D , where for each t ∈ {1, 2, . . . , T }, x t ∈ R D represents the tth observations (a.k.a., measurements) of all variables and x d t denotes the measurement of dth variable of x t . Let s t ∈ R denote the time-stamp when the tth observation is obtained and we assume that the first observation is made at 𝑴 = 1 1 0 1 0 1 1 0 1 1 0 0 0 1 𝚫 = 0.0 0.1 0.5 1.5 0.6 0.9 0.6 0.0 0.1 0.5 1.0 1.6 1.9 2.5 Figure 2 : An example of measurement vectors x t , time stamps s t , masking m t , and time interval δ t . time-stamp 0 (i.e., s 1 = 0). A time series X could have missing values. We introduce a masking vector m t ∈ {0, 1} D to denote which variables are missing at time step t. The masking vector for x t is given by

m d t = 1, if x d t is observed 0, otherwise

For each variable d, we also maintain the time interval δ d t ∈ R since its last observation as

δ d t =    s t -s t-1 + δ d t-1 , t > 1, m d t-1 = 0 s t -s t-1 , t > 1, m d t-1 = 1 0, t = 1

An example of these notations is illustrated in Figure 2 . In this paper, we are interested in the time series classification problem, where we predict the labels l n given the time series data D, where

D = {(X n , s n , M n , ∆ n )} N n=1

, and

EQUATION

Tn , s n = s

EQUATION

Tn , M n = m (n) 1 , . . . , m (n) Tn , ∆ n = δ (n) 1 , . . . , δ (n)

Tn , and l n ∈ {1, . . . , L}.

2.1. GRU-RNN FOR TIME SERIES CLASSIFICATION:

We investigate the use of recurrent neural networks (RNN) for time-series classification, as their recursive formulation allow them to handle variable-length sequences naturally. Moreover, RNN shares the same parameters across all time steps which greatly reduces the total number of parameters we need to learn. Among different variants of the RNN, we specifically consider an RNN with gated recurrent units (Cho et al., 2014; Chung et al., 2014) , but similar discussion and convolutions are also valid for other RNN models such as LSTM (Hochreiter & Schmidhuber, 1997) .

The structure of GRU is shown in Figure 3(a) . For each jth hidden unit, GRU has a reset gate r j t and an update gate z j t to control the hidden state h j t at each time t. The update functions are shown as follows:

r t = σ (W r x t + U r h t-1 + b r ) z t = σ (W z x t + U z h t-1 + b z ) ht = tanh (W x t + U (r t h t-1 ) + b) h t = (1 -z t ) h t-1 + z t ht

where matrices W z , W r , W , U z , U r , U and vectors b z , b r , b are model parameters. We use σ for element-wise sigmoid function, and for element-wise multiplication. This formulation assumes that all the variables are observed. A sigmoid or soft-max layer is then applied on the output of the GRU layer at the last time step for classification task.

Existing work on handling missing values lead to three possible solutions with no modification on GRU network structure. One straightforward approach is simply replacing each missing observation with the mean of the variable across the training examples. In the context of GRU, we have

x d t ← m d t x d t + (1 -m d t )x d (1) where xd = N n=1 Tn t=1 m d t,n x d t,n N n=1

Tn t=1 m d t,n . We refer to this approach as GRU-mean.

A second approach is exploiting the temporal structure in time series. For example, we may assume any missing value is same as its last measurement and use forward imputation (GRU-forward), i.e.,

x d t ← m d t x d t + (1 -m d t )x d t (2)

where t < t is the last time the d-th variable was observed.

Instead of explicitly imputing missing values, the third approach simply indicates which variables are missing and how long they have been missing as a part of input, by concatenating the measurement, masking and time interval vectors as

EQUATION

where x

(n) t can be either from Equation ( 1) or ( 2). We later refer to this approach as GRU-simple.

Several recent works (Lipton et al., 2016; Choi et al., 2015; Pham et al., 2016) use RNNs on EHR data to model diseases and to predict patient diagnosis from health care time series data with irregular time stamps or missing values, but none of them have explicitly attempted to capture and utilize the missing patterns into their RNNs via systematically modified network architectures. Choi et al. (2015) feeds medical codes along with its time stamps into GRU model to predict the next medical event. This feeding time stamps idea is equivalent to the baseline GRU-simple without feeding the masking, which we denote as GRU-simple (interval only). Pham et al. (2016) takes time stamps into LSTM model, and modify its forgetting gate by either time decay and parametric time both from time stamps. However, their non-trainable decay is not that flexible, and the parametric time also does not change RNN model structure and is similar to GRU-simple (interval only). In addition, neither of them consider missing values in time series medical records, and the time stamp input used in these two models is vector for one patient, but not matrix for each input variable of one patient as ours. Lipton et al. (2016) achieves their best performance on diagnosis prediction by feeding masking with zero-filled missing values. Their model is equivalent to GRU-simple without feeding the time interval, and no model structure modification is made for further capturing and utilizing missingness.

We denote their best model as GRU-simple (masking only). Conclusively, our GRU-simple baseline can be considered as a generalization from all related RNN models mentioned above and as shown in the experiments these GRU-simple variations have quite close performance.

These approaches solve the missing value issue to a certain extent, However, it is known that imputing the missing value with mean or forward imputation cannot distinguish whether missing values are imputed or truly observed. Simply concatenating masking and time interval vectors fails to exploit the temporal structure of missing values. Thus none of them fully utilize missingness in data to achieve desirable performance.

2.2. GRU-D: MODEL WITH TRAINABLE DECAYS:

To fundamentally address the issue of missing values in time series, we notice two important properties of the missing values in time series, especially in health care domains: First, the value of the missing variable tend to be close to some default value if its last observation happens a long time ago. This property usually exists in health care data for human body as homeostasis mechanisms and is considered to be critical for disease diagnosis and treatment (Vodovotz et al., 2013) . Second, the influence of the input variables will fade away over time if the variable has been missing for a while. For example, one medical feature in electronic health records (EHRs) is only significant in a certain temporal context (Zhou & Hripcsak, 2007) . Therefore we propose a GRU-based model called GRU-D, in which a decay mechanism is designed for the input variables and the hidden states to capture the aforementioned properties. We introduce decay rates in the model to control the decay mechanism by considering the following important factors. First, each input variable in health care time series has its own medical meaning and importance. The decay rates should be flexible to differ from variable to variable based on the underlying properties associated with the variables.

Second, as we see lots of missing patterns are informative in prediction tasks, the decay rate should be indicative of such patterns and benefits the prediction tasks. Furthermore, since the missing patterns are unknown and possibly complex, we aim at learning decay rates from the training data rather than being fixed a priori. That is, we model a vector of decay rates γ as

γ t = exp {-max (0, W γ δ t + b γ )} (4)

where W γ and b γ are model parameters that we train jointly with all the other parameters of the GRU. We chose the exponentiated negative rectifier in order to keep each decay rate monotonically decreasing in a reasonable range between 0 and 1. Note that other formulations such as a sigmoid function can be used instead, as long as the resulting decay is monotonic and is in the same range.

Our proposed GRU-D model incorporates two different trainable decays to utilize the missingness directly with the input feature values and implicitly in the RNN states. First, for a missing variable, we use an input decay γ x to decay it over time toward the empirical mean (which we take as a default configuration), instead of using the last observation as it is. Under this assumption, the trainable decay scheme can be readily applied to the measurement vector by

x d t ← m d t x d t + (1 -m d t )γ x d t x d t + (1 -m d t )(1 -γ x d t )x d (5)

where x d t is the last observation of the d-th variable (t < t) and xd is the empirical mean of the dth variable. When decaying the input variable directly, we constrain W γx to be diagonal, which effectively makes the decay rate of each variable independent from the others. Sometimes the input decay may not fully capture the missing patterns since not all missingness information can be represented in decayed input values. In order to capture richer knowledge from missingness, we also have a hidden state decay γ h in GRU-D. Intuitively, this has an effect of decaying the extracted features (GRU hidden states) rather than raw input variables directly. This is implemented by decaying the previous hidden state h t-1 before computing the new hidden state h t as

EQUATION

in which case we do not constrain W γ h to be diagonal. In addition, we feed the masking vectors (m t ) directly into the model. The update functions of GRU-D are

z t = σ (W z x t + U z h t-1 + V z m t + b z ) r t = σ (W r x t + U r h t-1 + V r m t + b r ) ht = tanh (W x t + U (r t h t-1 ) + V m t + b) h t = (1 -z t ) h t-1 + z t ht

where x t and h t-1 are respectively updated by Equation ( 5) and ( 6), and V z , V r , V are new parameters for masking vector m t .

To validate GRU-D model and demonstrate how it utilizes informative missing patterns, in Figure 4 , we show the input decay (γ x ) plots and hidden decay (γ h ) histograms for all the variables for predicting mortality on PhysioNet dataset. For input decay, we notice that the decay rate is almost constant for the majority of variables. However, a few variables have large decay which means that the model relies less on the previous observations for prediction. For example, the changes in the variable values of weight, arterial pH, temperature, and respiration rate are known to impact the ICU patients health condition. The hidden decay histograms show the distribution of decay parameters related to each variable. We noticed that the parameters related to variables with smaller missing rate are more spread out. This indicates that the missingness of those variables has more impact on decaying or keeping the hidden states of the models.

Notice that the decay term can be generalized to LSTM straightforwardly. In practical applications, missing values in time series may contain useful information in a variety of ways. A better model should have the flexibility to capture different missing patterns. In order to demonstrate the capacity of our GRU-D model, we discuss some model variations in Appendix A.2.

3. EXPERIMENTS:



3.1. DATASET DESCRIPTIONS AND EXPERIMENTAL DESIGN:

We demonstrate the performance of our proposed models on one synthetic and two real-world health-care datasets1 and compare it to several strong machine learning and deep learning approaches PhysioNet Challenge 2012 dataset (PhysioNet) This dataset, from PhysioNet Challenge 2012 (Silva et al., 2012) , is a publicly available collection of multivariate clinical time series from 8000 intensive care unit (ICU) records. Each record is a multivariate time series of roughly 48 hours and contains 33 variables such as Albumin, heart-rate, glucose etc. We used Training Set A subset in our experiments since outcomes (such as in-hospital mortality labels) are publicly available only for this subset. We conduct the following two prediction tasks on this dataset: 1) Mortality task: Predict whether the patient dies in the hospital. There are 554 patients with positive mortality label. We treat this as a binary classification problem. and 2) All 4 tasks: Predict 4 tasks: in-hospital mortality, length-of-stay less than 3 days, whether the patient had a cardiac condition, and whether the patient was recovering from surgery. We treat this as a multi-task classification problem.

MIMIC-III dataset (MIMIC-III):

This public dataset (Johnson et al., 2016) has deidentified clinical care data collected at Beth Israel Deaconess Medical Center from 2001 to 2012. It contains over 58,000 hospital admission records. We extracted 99 time series features from 19714 admission records for 4 modalities including input-events (fluids into patient, e.g., insulin), output-events (fluids out of the patient, e.g., urine), lab-events (lab test results, e.g., pH values) and prescription-events (drugs prescribed by doctors, e.g., aspirin). These modalities are known to be extremely useful for monitoring ICU patients. We only use the first 48 hours data after admission from each time series. We perform following two predictive tasks: 1) Mortality task: Predict whether the patient dies in the hospital after 48 hours. There are 1716 patients with positive mortality label and we perform binary classification. and 2) ICD-9 Code tasks: Predict 20 ICD-9 diagnosis categories (e.g., respiratory system diagnosis) for each admission. We treat this as a multi-task classification problem.

3.2. METHODS AND IMPLEMENTATION DETAILS:

We categorize all evaluated prediction models into three following groups:

• Non-RNN Baselines (Non-RNN): We evaluate logistic regression (LR), support vector machines (SVM) and Random Forest (RF) which are widely used in health care applications. • RNN Baselines (RNN): We take GRU-mean, GRU-forward, GRU-simple, and LSTM-mean (LSTM model with mean-imputation on the missing measurements) as RNN baselines. • Proposed Methods (Proposed): This is our proposed GRU-D model from Section 2.2. The non-RNN baselines cannot handle missing data directly. We carefully design experiments for non-RNN models to capture the informative missingness as much as possible to have fair comparison with the RNN methods. Since non-RNN models only work with fixed length inputs, we regularly sample the time-series data to get a fixed length input and perform imputation to fill in the missing values. Similar to RNN baselines, we can concatenate the masking vector along with the measurements and feed it to non-RNN models. For PhysioNet dataset, we sample the time series on an hourly basis and propagate measurements forward (or backward) in time to fill gaps. For MIMIC-III dataset, we consider two hourly samples (in the first 48 hours) and do forward (or backward) imputation.

Our preliminary experiments showed 2-hourly samples obtains better performance than one-hourly samples for MIMIC-III. We report results for both concatenation of input and masking vectors (i.e., SVM/LR/RF-simple) and only input vector without masking (i.e., SVM/LR/RF-forward). We use the scikit-learn (Pedregosa et al., 2011) for the non-RNN model implementation and tune the parameters by cross-validation. We choose RBF kernel for SVM since it performs better than other kernels.

For RNN models, we use a one layer RNN to model the sequence, and then apply a soft-max regressor on top of the last hidden state h T to do classification. We use 100 and 64 hidden units in GRU-mean for MIMIC-III and PhysioNet datasets, respectively. All the other RNN models were constructed to have a comparable number of parameters.2 For GRU-simple, we use mean imputation for input as shown in Equation (1). Batch normalization (Ioffe & Szegedy, 2015) and dropout (Srivastava et al., 2014) of rate 0.5 are applied to the top regressor layer. We train all the RNN models with the Adam optimization method (Kingma & Ba, 2014) and use early stopping to find the best weights on the validation dataset. All the input variables are normalized to be 0 mean and 1 standard deviation. We report the results from 5-fold cross validation in terms of area under the ROC curve (AUC score). We provide more detailed comparisons of RNN baselines and variations in Appendix A.3.4 and evaluations on multilayer RNN models in Appendix A.3.5.

3.3. QUANTITATIVE RESULTS:

Exploiting informative missingness on synthetic dataset As illustrated in Figure 1 , missing patterns can be useful in solving prediction tasks. A robust model should exploit informative missingness properly and avoid inducing nonexistent relations between missingness and predictions. To evaluate the impact of modeling missingness we conduct experiments on the synthetic Gesture datasets. We process the data in 4 different settings with the same missing rate but different correlations between missing rate and the label. A higher correlation implies more informative missingness. Figure 5 shows the AUC score comparison of three GRU baseline models (GRU-mean, GRU-forward, GRU-simple) and the proposed GRU-D. Since GRU-mean and GRU-forward do not utilize any missingness (i.e., masking or time interval), they perform similarly across all 4 settings. GRU-simple and GRU-D benefit from utilizing the missingness, especially when the correlation is high. Our GRU-D achieves the best performance in all settings, while GRU-simple fails when the correlation is low. The results on synthetic datasets demonstrates that our proposed model can model and distinguish useful missing patterns in data properly compared with baselines.

Prediction task evaluation on real datasets:

We evaluate all methods in Section 3.2 on MIMIC-III and PhysioNet datasets. We noticed that dropout in the recurrent layer helps a lot for all RNN models

Under review as a conference paper at ICLR 2017 on both of the datasets, probably because they contain more input variables and training samples than synthetic dataset. Similar to Gal (2015) , we apply dropout rate of 0.3 with same dropout samples at each time step on weights W , U , V . Table 2 shows the prediction performance of all the models on mortality task. All models except for random forest improve their performance when they feed missingness indicators along with inputs. The proposed GRU-D achieves the best AUC score on both datasets. We also conduct multi-task classification experiments for all 4 tasks on PhysioNet and 20 ICD-9 code tasks on MIMIC-III using all the GRU models. As shown in Table 1 , GRU-D performs best in terms of average AUC score across all tasks and in most of the single tasks. 

3.4. DISCUSSIONS:

Online prediction in early stage Although our model is trained on the first 48 hours data and makes prediction at the last time step, it can be used directly to make predictions before it sees all the time series and can make predictions on the fly. This is very useful in applications such as health care, where early decision making is beneficial and critical for patient care. Figure 6 shows the online prediction results for MIMIC-III mortality task. As we can see, AUC is around 0.7 at first 12 hours for all the GRU models and it keeps increasing when longer time series is fed into these models. GRU-D and GRU-simple, which explicitly handle missingness, perform consistently better than the other two methods. In addition, GRU-D outperforms GRU-simple when making predictions given time series of more than 24 hours, and has at least 2.5% higher AUC score after 30 hours. This indicates that GRU-D is able to capture and utilize long-range temporal missing patterns. Furthermore, GRU-D achieves similar prediction performance (i.e., same AUC) as best non-RNN baseline model with less time series data. As shown in the figure, GRU-D has same AUC performance at 36 hours as the best non-RNN baseline model (RF-simple) at 48 hours. This 12 hour improvement of GRU-D over non-RNN baseline is highly significant in hospital settings such as ICU where time-saving critical decisions demands accurate early predictions.

Model Scalability with growing data size In many practical applications, model scalability with large dataset size is very important. To evaluate the model performance with different training dataset size, we subsample three smaller datasets of 2000 and 10000 admissions from the entire MIMIC-III dataset while keeping the same mortality rate. We compare our proposed models with all GRU baselines and two most competitive non-RNN baselines (SVM-simple, RF-simple) and shows the prediction results in Figure 7 . We observe that all models can achieve improved performance given more training samples. However, the improvements of non-RNN baselines are quite limited compared to GRU models, and our GRU-D model achieves the best results on the two larger datasets. These results indicate the performance gap between GRU-D and non-RNN baselines will continue to grow as more data become available.

4. SUMMARY:

In this paper, we proposed novel GRU-based model to effectively handle missing values in multivariate time series data. Our model captures the informative missingness by incorporating masking and time interval directly inside the GRU architecture. Empirical experiments on both synthetic and real-world health care datasets showed promising results and provided insightful findings. In our future work, we will explore deep learning approaches to characterize missing-not-at-random data and we will conduct theoretical analysis to understand the behaviors of existing solutions for missing values.

A APPENDIX A.1 INVESTIGATION OF RELATION BETWEEN MISSINGNESS AND LABELS:

In many time series applications, the pattern of missing variables in the time series is often informative and useful for prediction tasks. Here, we empirically confirm this claim on real health care dataset by investigating the correlation between the missingness and prediction labels (mortality and ICD-9 diagnosis categories). We denote the missing rate for a variable d as p d X and calculate it by

p d X = 1 -1 T T t=1 m d t .

Note that p d X is dependent on mask vector (m d t ) and number of time steps T . For each prediction task, we compute the Pearson correlation coefficient between p d X and label across all the time series. As shown in Figure 1 , we observe that on MIMIC-III dataset the missing rates with low rate values are usually highly (either positive or negative) correlated with the labels. The distinct correlation between missingness and labels demonstrates usefulness of missingness patterns in solving prediction tasks.

A.2 GRU-D MODEL VARIATIONS:

In this section, we will discuss some variations of GRU-D model, and also compare some related RNN models which are used for time series with missing data with the proposed model. 

A.2.1 GRU MODEL WITH DIFFERENT TRAINABLE DECAYS:

The proposed GRU-D applies trainable decays on both input and hidden state transitions in order to capture the temporal missing patterns explicitly. This decay idea can be straightforwardly generated to other parts inside the GRU models separately or jointly, given different assumptions on the impact of missingness. As comparisons, we also describe and evaluate several modifications of GRU-D model. 8 (a)) and GRU-DS (Figure 8(b) ) decay only the input and only the hidden state by Equation ( 5) and ( 6), respectively. They can be considered as two simplified models of the proposed GRU-D. GRU-DI aims at capturing direct impact of missing values in the data, while GRU-DS captures more indirect impact of missingness. Another intuition comes from this perspective: if an input variable is just missing, we should pay more attention to this missingness; however, if an variable has been missing for a long time and keeps missing, the missingness becomes less important. We can utilize this assumption by decaying the masking. This brings us the model GRU-DM shown in Figure 8(c) , where we replace the masking m d t fed into GRU-D in by

GRU-DI (Figure:

EQUATION

where the equality holds since m d t is either 0 or 1. We decay the masking for each variable independently from others by constraining W γm to be diagonal.

A.2.2 GRU-IMP: GOAL-ORIENTED IMPUTATION MODEL:

We may alternatively let the GRU-RNN predict the missing values in the next timestep on its own. When missing values occur only during test time, we simply train the model to predict the measurement vector of the next time step as a language model (Mikolov et al., 2010) and use it to fill the missing values during test time. This is unfortunately not applicable for some time series applications such as in health care domains, which also have missing data during training.

Instead, we propose goal-oriented imputation model here called GRU-IMP, and view missing values as latent variables in a probabilistic graphical model. Given a timeseries X, we denote all the missing variables by M X and all the observed ones by O X . Then, training a time-series classifier with missing variables becomes equivalent to maximizing the marginalized log-conditional probability of a correct label l, i.e., log p(l|O X ).

The exact marginalized log-conditional probability is however intractable to compute, and we instead maximize its lowerbound:

log p(l|O X ) = log M X p (l|M X , O X ) p (M X |O X ) ≥ E M X ∼p(M X |O X ) log p (l|M X , O X )

where we assume the distribution over the missing variables at each time step is only conditioned on all the previous observations:

p (M X |O X ) = T t=1 m d t =1

1≤d≤D p(x d t |x 1:(t-1) , m 1:(t-1) , δ 1:(t-1) )

Although this lowerbound is still intractable to compute exactly, we can approximate it by Monte Carlo method, which amounts to sampling the missing variables at each time as the RNN reads the input sequence from the beginning to the end, such that

x d t ← m d t x d t + (1 -m d t )x d t (9)

where xt ∼ x d t |x 1:(t-1) , m 1:(t-1) , δ 1:(t-1) . By further assuming that xt ∼ N µ t , σ 2 t , µ t = γ t (W x h t-1 + b x ) and σ t = 1, we can use a reparametrization technique widely used in stochastic variational inference (Kingma & Welling, 2013; Rezende et al., 2014) to estimate the gradient of the lowerbound efficiently. During the test time, we simply use the mean of the missing variable, i.e., xt = µ t , as we have not seen any improvement from Monte Carlo approximation in our preliminary experiments. We view this approach as a goal-oriented imputation method and show its structure in Figure 8(d ). The whole model is trained to minimize the classification cross-entropy error log loss and we take the negative log likelihood of the observed values as a regularizer. For each of the three datasets used in our experiments, we list the number of samples, the number of input variables, the mean and max number of time steps for all the samples, and the mean of all the variable missing rates in Table 3 . 

:

A summary statistics of the three datasets is shown in Appendix A.3.1.

Appendix A.3.2 compares all GRU models tested in the experiments in terms of model size.

Figures and Tables:

Figure 1: Demonstrations of informative missingness on MIMIC-III dataset. Left figure shows variable missing rate (x-axis, missing rate; y-axis, input variable). Middle/right figures respectively shows the correlations between missing rate and mortality/ICD-9 diagnosis categories (x-axis, target label; y-axis, input variable; color, correlation value). Please refer to Appendix A.1 for more details.

Figure 3: Graphical illustrations of the original GRU (left) and the proposed GRU-D (right) models.

Figure 4: Plots of input decay γ xt (top) and histrograms of hidden state decay γ ht (bottom) of all 33 variables in GRU-D model for predicting mortality on PhysioNet dataset. Variables in green are lab measurements; variables in red are vital signs; mr refers to missing rate.

Figure5: Classification performance on Gesture synthetic datasets. x-axis: average Pearson correlation of variable missing rates and target label in that dataset; y-axis: AUC score.

Figure 8: Graphical illustrations of variations of proposed GRU models.

-axis, value of decay parameter W γ h ; y-axis, count.

Table content: <table><tr><td>10 1 10 2</td><td colspan="2">32: pH mr: 0.9118</td><td colspan="2">9: FiO2 mr: 0.883</td><td colspan="2">23: RespRate mr: 0.8053</td><td colspan="2">10: GCS mr: 0.7767</td><td colspan="2">26: Temp mr: 0.6915</td><td colspan="2">31: Weight mr: 0.5452</td><td colspan="2">29: Urine mr: 0.5095</td><td colspan="2">17: MAP mr: 0.2141</td><td colspan="2">8: DiasABP mr: 0.2054</td><td colspan="2">25: SysABP mr: 0.2052</td><td colspan="2">14: HR mr: 0.1984</td></tr><tr><td/><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td></tr><tr><td/><td/><td/><td/><td/><td colspan="2">(b) x</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

Model performances measured by average AUC score (mean ± std) for multi-task predictions on real datasets. Results on each class are shown in Appendix A.3.3 for reference.

Table content: <table><tr><td>Models</td><td>MIMIC-III ICD-9 20 tasks</td><td>PhysioNet All 4 tasks</td></tr><tr><td>GRU-mean</td><td>0.7070 ± 0.001</td><td>0.8099 ± 0.011</td></tr><tr><td>GRU-forward</td><td>0.7077 ± 0.001</td><td>0.8091 ± 0.008</td></tr><tr><td>GRU-simple</td><td>0.7105 ± 0.001</td><td>0.8249 ± 0.010</td></tr><tr><td>GRU-D</td><td colspan="2">0.7123 ± 0.003 0.8370 ± 0.012</td></tr></table>

Model performances measured by AUC score (mean ± std) for mortality prediction.

Table content: <table><tr><td colspan="2">Models</td><td>MIMIC-III</td><td>PhysioNet</td></tr><tr><td/><td>LR-forward</td><td>0.7589 ± 0.015</td><td>0.7423 ± 0.011</td></tr><tr><td/><td>SVM-forward</td><td>0.7908 ± 0.006</td><td>0.8131 ± 0.018</td></tr><tr><td>Non-RNN</td><td>RF-forward</td><td>0.8293 ± 0.004</td><td>0.8183 ± 0.015</td></tr><tr><td/><td>LR-simple</td><td>0.7715 ± 0.015</td><td>0.7625 ± 0.004</td></tr><tr><td/><td>SVM-simple</td><td>0.8146 ± 0.008</td><td>0.8277 ± 0.012</td></tr><tr><td/><td>RF-simple</td><td>0.8294 ± 0.007</td><td>0.8157 ± 0.013</td></tr><tr><td/><td>LSTM-mean</td><td>0.8142 ± 0.014</td><td>0.8025 ± 0.013</td></tr><tr><td>RNN</td><td>GRU-mean</td><td>0.8192 ± 0.013</td><td>0.8195 ± 0.004</td></tr><tr><td/><td>GRU-forward</td><td>0.8252 ± 0.011</td><td>0.8162 ± 0.014</td></tr><tr><td/><td>GRU-simple</td><td>0.8380 ± 0.008</td><td>0.8155 ± 0.004</td></tr><tr><td colspan="2">Proposed GRU-D</td><td colspan="2">0.8527 ± 0.003 0.8424 ± 0.012</td></tr></table>

Dataset statistics.

Table content: <table><tr><td/><td colspan="3">MIMIC-III PhysioNet2012 Gesture</td></tr><tr><td># of samples (N )</td><td>19714</td><td>4000</td><td>378</td></tr><tr><td># of variables (D)</td><td>99</td><td>33</td><td>23</td></tr><tr><td>Mean of # of time steps</td><td>35.89</td><td>68.91</td><td>21.42</td></tr><tr><td>Maximum of # of time steps</td><td>150</td><td>155</td><td>31</td></tr><tr><td>Mean of variable missing rate</td><td>0.9621</td><td>0.8225</td><td>N/A</td></tr></table>

Peer review (iclr2017):
This paper proposed a way to deal with supervised multivariate time series tasks involving missing values. The high level idea is still using the recurrent neural network (specifically, GRU in this paper) to do sequence supervised learning, e.g., classification, but modifications have been made to the input and hidden layers of RNNs to tackle the missing value problem. pros: 1) the insight of utilizing missing value is critical. the observation of decaying effect in the healthcare application is also interesting; 2) the experiment seems to be solid; the baseline algorithms and analysis of results are also done properly. cons: 1) the novelty of this work is not enough. Adding a decaying smooth factor to input and hidden layers seems to be the main modification of the architecture. 2) the datasets used in this paper are small. 3) the decaying effect might not be able to generalize to other domains.