GRADIENT IMPORTANCE LEARNING FOR INCOMPLETE OBSERVATIONS
Year: 2022

Abstract:
Though recent works have developed methods that can generate estimates (or imputations) of the missing entries in a dataset to facilitate downstream analysis, most depend on assumptions that may not align with real-world applications and could suffer from poor performance in subsequent tasks such as classification. This is particularly true if the data have large missingness rates or a small sample size. More importantly, the imputation error could be propagated into the prediction step that follows, which may constrain the capabilities of the prediction model. In this work, we introduce the gradient importance learning (GIL) method to train multilayer perceptrons (MLPs) and long short-term memories (LSTMs) to directly perform inference from inputs containing missing values without imputation. Specifically, we employ reinforcement learning (RL) to adjust the gradients used to train these models via back-propagation. This allows the model to exploit the underlying information behind missingness patterns. We test the approach on real-world time-series (i.e., MIMIC-III), tabular data obtained from an eye clinic, and a standard dataset (i.e., MNIST), where our imputation-free predictions outperform the traditional two-step imputation-based predictions using state-of-the-art imputation methods.

heading: 1 INTRODUCTION
text: We consider learning from incomplete datasets. Components of the data could be missing for various reasons, for instance, missing responses from participants, data loss, and restricted access issues. This phenomenon is prevalent in the healthcare domain, mostly in the context of electronic health records (EHRs), which are structured as patient-specific irregular timelines with attributes, e.g., diagnosis, laboratory tests, vitals, thus resulting in high missingness across patients for any arbitrary time point. Such missingness introduces difficulties when developing models and performing inference in real-world applications such as Kam & Kim (2017); Scherpf et al. (2019). Existing works tackle this problem by either proposing imputation algorithms (IAs) to explicitly produce estimates of the missing data, or by imposing imputation objectives during inference, e.g., withholding observed values as being the ground-truth and learning to impute them. However, some of these either require additional modeling assumptions for the underlying distributions (Bashir & Wei, 2018; Fortuin et al., 2020), or formatting of the data (Yu et al., 2016; Schnabel et al., 2016). Other applications depend on domain knowledge for pre-processing and modeling such as Yang et al. (2018); Kam & Kim (2017), or introduce additional information-based losses (Cao et al., 2018; Lipton et al., 2016), which are usually intractable with real-world data. Moreover, generative methods (Mattei & Frellsen, 2019; Yoon et al., 2018) could result in imputations with high variation (e.g., low confidence of the output distribution), when data have high missingness rates or small sample sizes. Figure 1 illustrates imputations generated by two state-of-the-art deep generative models on the MNIST dataset with 50%, 70% and 90% of the pixels set as missing (i.e., masked out). It is observed that both approaches suffer from inaccurate reconstruction of the original digits as the missingness rate increases, which is manifested as some of the imputed images not being recognizable as digits. Importantly, the error introduced by the imputation process can be further propagated into downstream inference and prediction stages.
Imputing the missing data (either explicitly or implicitly) is, in most cases, not necessary, more so, considering that sometimes where and when the data are missing can be intrinsically informative. ∗Duke University, USA. Contact: {qitong.gao, miroslav.pajic}@duke.edu. †King Abdullah University of Science and Technology, Saudi Arabia. Code available at https://github.com/gaoqitong/gradient-importance-learning.
Consider a scenario in which two patients, A and B, admitted to the intensive care unit (ICU) suffer from bacterial and viral infections, respectively, and assume that the healthcare provider monitors the status of the patients by ordering two (slightly) different blood tests periodically, namely, a culture test/panel specific to bacterial infections and a RT-PCR test/panel specific to viral infections. Hence, patient A is likely to have much fewer orders (and results) for viral tests. In both cases, the missingness patterns are indicative of the underlying condition of both patients. Moreover, such patterns are more commonly found in incomplete data caused by missing at random (MAR) or missing not at random (MNAR) mechanisms, which usually introduce additional difficulties for inference (Little & Rubin, 2019). Inspired by this, we propose the gradient importance learning (GIL) method, which facilitates an imputation-free learning framework for incomplete data by simultaneously leveraging the observed data and the information underlying missingness patterns.
Specifically, we propose to use an importance matrix to weight the gradients that are used to update the model parameters in the backpropagation process, after the computational graph of the model is settled, which allows the model to exploit the information underlying missingness without imputation. However, these gradients cannot be tracked by automated differentiation tools such as Tensorflow (Abadi et al., 2016). So motivated, we propose to generate the importance matrix using a reinforcement learning (RL) policy on-the-fly. This is done by conditioning on the status of training procedure
characterized by the model parameters, inputs and model performance at the current training step. Concurrently, RL algorithms are used to update the policy by first modeling the back-propagation process as an RL environment, or a Markov decision process (MDP). Then, by interacting with the environment, the RL algorithm can learn the optimal policy so the importance matrix can aid the training of the prediction models to attain desirable performance. Moreover, we also show that our framework can be augmented with feature learning techniques, e.g., contrastive learning as in Chen et al. (2020), which can further improve the performance of the proposed imputation-free models.
The technical contributions of this work can be summarized as follows: (i) to the best of our knowledge, this is the first work that trains deep learning (DL) models to perform accurate imputationfree predictions with missing inputs. This allows models to effectively handle high missing rates, while significantly reducing the prediction error compared to existing imputation-prediction frameworks. (ii) Unlike existing approaches that require additional modeling assumptions or that rely on methods (or pre-existing knowledge) that intrinsically have advantages over specific types of data, our method does not require any assumptions or domain expertise. Consequently, it can consistently achieve top performance over a variety of data and under different conditions. (iii) The proposed framework also facilitates feature learning from incomplete data, as the importance matrix can guide the hidden layers of the model to capture information underlying missingness patterns during training, which results in more expressive features, as illustrated in Figure 3.

heading: 2 GRADIENT IMPORTANCE LEARNING (GIL)
text: In this section, we introduce GIL; a method that facilitates training of imputation-free prediction models with incomplete data. In Section 2.2, we show how the gradient descent directions can be adjusted using the importance matrix A. In Section 2.3, we introduce an RL framework to generate the elements of A which can leverage the information underlying missingness during training, as illustrated in Figure 2.

heading: 2.1 PROBLEM FORMULATION
text: We consider a dataset containing N observations X = (X1,X2, . . . ,XN ) where each Xj can be represented by a vector xj ∈ Rd (for tabular data) or a matrix Xj ∈ RTj×d (for time-series) with
Tj denoting the time horizon of the sequence Xj . Since the focus of this work is principally on time-series data, recurrent neural networks will receive significant attention. We also define the set of missing indicatorsM = (M1,M2, . . . ,MN ), where Mj = mj ∈ {0, 1}d or Mj ∈ {0, 1}Tj×d depending on the dimension of Xj , and 0’s (or 1’s) correspond to the entries in Xj that are missing (or not missing), respectively. Finally, we assume that each Xj is associated with some label yj ∈ Y ; thus Y = (y1,y2, . . . ,yN ) denotes the set of labels. We define the problem as learning a model parameterized by W to directly predict yj by generating ŷj that maximizes the log-likelihood∑N j=1 log p(ŷj |Xj ,Mj ,W), without imputing the missing values in Xj .

heading: 2.2 GRADIENT IMPORTANCE
text: Missingness is a common issue found in tabular data and time series, where multi-layered perceptron (MLP) (Bishop, 2006) and long short-term memory (LSTM) (Hochreiter & Schmidhuber, 1997) models, respectively, are often considered for predictions. Below we illustrate how the importance matrix can be applied to gradients, which are produced by taking regular stochastic gradient descent (SGD) steps, toward training MLP models with incomplete tabular inputs. Note that this idea can be easily extended to sequential inputs with LSTM models and the corresponding details can be found in Appendix A. First, the definitions of the dataset and missing indicators can be reduced to Xtab = (x1,x2, . . . ,xN ) and Mtab = (m1,m2, . . . ,mN ), with xj ∈ Rd and mj ∈ Rd, respectively. Then we define an MLP with k hidden layers, for k ≥ 2, as follows
ŷ = φout(Woutφk(Wk . . . φ2(W2φ1(W1x)))), (1)
where x ∈ Rd is the input, Wi and φi are the weight matrix and activation functions for the i-th hidden layer, respectively, and we have omitted the bias terms for notational convenience. The first layer can be interpreted as the encoding layer Wenc = W1 that maps inputs to features ζ = fenc(x|Wenc), while the rest are the inference layers Winf = {W2, . . . ,Wk,Wout} that map the features to prediction ŷ = finf (ζ|Winf ). In the following proposition, we show that the gradients of some commonly used loss functions, E(ŷ,y), such as cross entropy or mean squared errors, w.r.t. W1, can be formulated in the form of an outer product. The proof and details can be found in Appendix E.
Proposition 1 Given a MLP and a smooth loss function E(ŷ,y), the gradients of E w.r.t. the encoding layer can be formulated as ∂E/∂W1 = ∆1x>, where ∆1 contains the gradients propagated from all the inference layers, x ∈ Xtab and y ∈ Y .
From this Proposition it can be seen that the gradients that are used to train the encoding layers, following regular SGD solvers, can be formulated as the outer product between the gradients ∆ propagated from the inference layers and the input x as
(∂E/∂Wenc)SGD = ∆ · x>, (2)
where Wenc ∈ Re×d, ∆ ∈ Re, x ∈ Rd, e is the dimension of the features, and d is the dimension of the input. To simplify notation, note that henceforth we use x to refer to an individual tabular data xj ∈ Xtab, regardless of its index. The corresponding SGD updates for training Wenc, using learning rate α, are Wenc ←Wenc − α∆ · x>. As a result, it is observed from (2) that the j-th (j ∈ [1, d]) column of (∂E/∂Wenc)SGD is weighted by the j-th element in x given ∆. However, given that some entries in x could be missing (according to the missing indicator m) and their values are usually replaced by a placeholder ξ ∈ R, it may not be ideal do directly back-propagate the gradients in the form of (2), as the features ζ captured by the encoding layers may not be expressive enough for the inference layers to make accurate predictions. Instead, we consider introducing the importance matrix A ∈ [0, 1]e×d to adjust the gradient descent direction as
Wenc ←Wenc − α(∂E/∂Wenc)SGD A, (3) which not only trains the model to capture information from the observed inputs that are most useful for making accurate predictions, but also to leverage the information underlying the missingness in the data. The elements of A can be generated using the RL framework introduced in the next section.
Note that above we have omitted all the bias terms to simplify our presentation and note that: i) gradients w.r.t. to biases do not conform to the outer product format, and ii) more importantly, these gradients do not depend on the inputs – thus in practice, the importance matrix A is only applied
to the gradients of the weights. Moreover, we do not consider convolutional neural network (CNN) models (LeCun et al., 1999) in this work because of i). Though the proposed framework could still be applied to CNNs, it may not be as efficient as for MLPs or LSTMs, where the search space for A is significantly reduced by taking advantages of the outer product format, as shown in (4) below. Importantly, our focus is to address the missingness in tabular data and time-series, in which case MLPs and LSTMs are appropriate, while the missingness in images is usually related to applications in other domains such as compressed sensing (Wu et al., 2019) or image super-resolution (Dong et al., 2015), which we plan to explore in the future.

heading: 2.3 RL TO GENERATE IMPORTANCE MATRIX
text: Now we show how to use RL to generate the importance matrix A. In general, RL aims to learn an optimal policy that maximizes expected rewards, by interacting with an unknown environment (Sutton & Barto, 2018). Formally, the RL environment is characterized by a Markov decision process (MDP), with details introduced below. Each time after the agent chooses an action a = π(s) at state s following some policy π, the environment responds with the next state s′, along with an immediate reward r = R(s, a) given a reward function R(·, ·). The goal is to learn an optimal policy that maximizes the expected total reward defined as J(π) = Eρπ [ ∑T i=0 γ
ir(si, ai)] , where T is the horizon of the MDP, γ ∈ [0, 1) is the discounting factor, and ρπ = {(s0, a0), (s1, a1), . . . |ai = π(si)} is the sequence of states and actions drawn from the trajectory distribution determined by π. In our case, the elements of A are determined on-the-fly following an RL policy πθ, parameterized by θ, conditioned on some states that characterize the status of the back-propagation at the current time step. The policy πθ is updated, concurrently with the weights Wenc, by an RL agent that interacts with the back-propagation process modeled as the MDP defined asM = (S,A,P, R, γ), with each of its elements introduced below.
State Space S. The state space characterizes the agent’s knowledge of the environment at each time step. To constitute the states we consider the 4-tuple s = (x,m, ζ, ŷ) including current training input x ∈ Rd, the missing indicator m ∈ Rd, the feature ζ, i.e., the outputs of the encoding layer of MLPs or the hidden states ht of LSTMs, and the predictions ŷ produced by the inference layers.
Action Space A. We consider combining the importance matrix A ∈ [0, 1]e×d with the parameter gradients (∂E/∂Wenc)SGD when updating Wenc following W′enc ←Wenc − α(∂E/∂Wenc)SGD A = Wenc − α∆ · (x> a>); (4) where here the equality follows from (2) and Proposition 1, such that all rows of A can be set to be equal to the importance a> ∈ [0, 1]d, which is obtained from the policy following a = πθ(s). Transitions P : S × A → S . The transition dynamics determines how to transit from a current state s to the next state s′ given an action a in the environment. In our case, the pair (x, m) is sampled from the training dataset at each step, so the transitions x → x′ and m → m′ are determined by how training samples are selected from the training batch during back-propagation. The update of Wenc → W′enc is conditioned on the importance a as captured in (4), which results in the transitions ζ = fenc(x|Wenc) → ζ ′ = fenc(x′|W′enc). The update of Winf → W′inf follows from the regular SGD updates. Then, the transition ŷ = f(x|W) → ŷ′ = f(x′|W′) follows from the other transitions defined above. Finally, we can define the transition between states as s = (x,m, ζ, ŷ)→ s′ = (x′,m′, ζ ′, ŷ′). Reward Function R. After the RL agent takes an action a at the state s, an immediate reward r = R(s,a) is returned by the environment which provides essential information to update πθ (Silver et al., 2014; Lillicrap et al., 2015). We define the reward function as R(s,a) = −E(ŷ,y). We utilize actor-critic RL (Silver et al., 2014; Lillicrap et al., 2015) to update πθ, which outputs the importance a that is used to concurrently update Wenc. Specifically, we train the target policy (or actor) πθ along with the critic Qν : S × A → R, parameterized by ν, by maximizing Jβ(πθ) = Es∼ρβ [Qν(s, πθ(s))] which gives an approximation to the expected total reward J(π). Specifically, the trajectories ρβ = {(s0, a0), (s1, a1), . . . |ai = β(si)} collected under the behavioral policy β : S → A are used to update θ and ν jointly following ν′ ← ν + ανδ∇νQν(s,a), θ′ ← θ + αθ∇θπθ(s)∇aQν(s,a)|a=πθ(s); (5) where β is usually obtained by adding noise to the output of πθ to ensure sufficient exploration of the state and action space, δ = r+ γQν(s′, πθ(s′))−Qν(s,a) is the temporal difference error (residual) in RL, and αθ, αν are the learning rates for θ, ν, respectively.
Algorithm 1 Gradient Importance Learning (GIL). Input: X , Y ,M, Wenc, Winf , πθ, Qν , αθ, αν , α, E Begin:
1: Initialize Wenc and Winf , actor πθ and critic Qν 2: Sample x from X and obtain the corresponding label y from Y 3: Obtain the feature ζ ← fenc(x|Wenc) and prediction ŷ = finf (ζ|Winf ) from the encoding
and inference layers, respectively 4: s← (x,m, ζ, ŷ) 5: for iter in 1 : max_iter do 6: Obtain importance from a behavioral policy a = β(s|πθ) 7: Train the encoding layer following W′enc ←Wenc − α∆ · (x> a>) as in (4) 8: Train the inference layers following regular gradient descent, i.e., W′inf ←Winf − α(∂E/∂Winf )SGD 9: Obtain the prediction following the updated weights ŷ← f(x|W′enc,W′inf )
10: Obtain the reward r ← R(s,a) 11: Get a new sample x′ from X and obtain the corresponding label y′ from Y 12: Obtain the feature ζ ′ ← fenc(x′|W′enc) and prediction ŷ′ = finf (ζ ′|W′inf ) from the encoding and inference layers, respectively 13: s′ ← (x′,m′, ζ ′, ŷ′) 14: Update the actor πθ and critic Qν using (s,a, r, s′) following (5) 15: s← s′,Wenc ←W′enc,Winf ←W′inf 16: end for
We summarize the GIL approach in Algorithm 1 and the detailed descriptions can be found in Appendix B. Note that the missing indicator m ∈ Rd would be a good heuristic to replace the importance a as it could prevent the gradients produced by the missing dimensions from propagation. However, it does not train the model to capture the hidden information underlying the missing data, which results in its performance being dominated by GIL as demonstrated in Section 4.

heading: 2.4 EXTENSIONS OF THE FRAMEWORK
text: The proposed GIL framework uses RL agents to guide the model toward minimizing its objective E(ŷ,y), which is usually captured by general prediction losses such as cross entropy or mean squared error. Below we use an example to show how our method can be extended to adapt ideas from related domains, which can augment training of the imputation-free prediction models by altering the reward function. Recent works in contrastive learning, e.g., Chen et al. (2020), can train DL models to generate highly expressive features by penalizing (or promoting) the distributional difference D(ζ+, ζ−) among the features associated with inputs that are similar to (or different from) each other, where ζ+ denotes the set of features corresponding to one set of homogeneous inputs Xcon, and ζ− are the features generated by data that are significantly different from the ones in Xcon. However, such methods may not be directly applied to the missing data setting considered in this work, as their loss D is usually designed toward unsupervised training and might need to be carefully re-worked in our case. This idea can be adapted to our framework by defining ζ+ as the features mapped from the inputs (by the encoding layers Wenc) associated with the same label y ∈ Y and ζ− the features corresponding to some other label y′ ∈ Y . Then we can define the new reward function R(s,a) = −E(ŷ,y) + c · D(ζ+, ζ−), c > 0, which does not require D to be carefully crafted, provided that Wenc is not trained by directly propagating gradients from it. In Section 4 it will be shown in case studies that by simply defining D as the mean squared error between ζ+ and ζ− can improve the prediction performance.

heading: 3 RELATED WORK
text: Missing Data Imputation. Traditional mean/median-filling and carrying-forward imputation methods are still used widely, as they are straightforward to implement and interpret (Honaker & King, 2010). Recently, there have been state-of-the-art imputation algorithms (IAs) proposed to produce smooth imputations with interpretable uncertainty estimates. Specifically, some adopt Bayesian methods where the observed data are fit to data-generating models, including Gaussian processes (Wilson et al., 2016; Fortuin & Rätsch, 2019; Fortuin et al., 2018), multivariate imputation by chained equations (MICE) (Azur et al., 2011), random forests (Stekhoven & Bühlmann, 2012), etc., or statistical
optimization methods such as expectation-maximization (EM) (García-Laencina et al., 2010; Bashir & Wei, 2018). However, they often suffer from limited scalability, require assumptions over the underlying distribution, or fail to generalize well when used with datasets containing mixed types of variables, i.e., when discrete and continuous values exist simultaneously (Yoon et al., 2018; Fortuin et al., 2020). There also exist matrix-completion methods that usually assume the data are static, i.e., information does not change over time, and often rely on low-rank assumptions (Wang et al., 2006; Yu et al., 2016; Mazumder et al., 2010; Schnabel et al., 2016). More recently, several DL-based IAs have been proposed following advancements in deep generative models such as deep latent variable models (DLVMs) and generative adversarial networks (GANs) (Mattei & Frellsen, 2019; Yoon et al., 2018; Fortuin et al., 2020; Li et al., 2019; Ma et al., 2018). DLVMs are usually trained to capture the underlying distribution of the data, by maximizing an evidence lower bound (ELBO), which could introduce high variations to the output distributions (Kingma & Welling, 2013) and cause poor performance in practice, if the data have high missingness (as for the example shown in Figure 1). There also exist end-to-end methods that withhold observed values from inputs and impose reconstruction losses during training of prediction models (Cao et al., 2018; Ipsen et al., 2020). In addition, data-driven IAs are developed specifically toward medical time series by incorporating prior domain knowledge (Yang et al., 2018; Scherpf et al., 2019; Calvert et al., 2016; Gao et al., 2021). On the other hand, a few recent works attempt to address missing inputs through an imputation-free manner (Morvan et al., 2020; Sportisse et al., 2020), however, they are limited to linear regression.
Attention. The importance matrix used in the GIL is somewhat reminiscent of visual attention mechanisms in the context of image captioning or multi-object recognition (Mnih et al., 2014; Ba et al., 2014; Xu et al., 2015) as they both require training of the prediction models with RL. We briefly discuss the connections and distinctions between them, with full details provided in Appendix C. Visual attentions are commonly used to train CNNs to focus on specific portions of the inputs that are most helpful for making predictions. They are determined by maximizing an evidence lower bound (ELBO), which is later proved to be equivalent to the REINFORCE algorithm in the RL literature (Sutton & Barto, 2018; Mnih et al., 2014). However, these methods cannot be applied directly to our problem, as they require features to be exclusively associated spatially with specific parts of the inputs, which is attainable by using convolutional encoders with image inputs but intractable with the type of data considered in our work. Instead, our approach overcomes this issue by directly applying importance weights, generated by RL, into the gradient space during back-propagation. Such issues motivate use of the term importance instead of attention. Lastly, our method does not require one to formulate the learning objective as an ELBO, and as a result GIL can adopt any state-of-the-art RL algorithm, without being limited to REINFORCE as in Mnih et al. (2014); Ba et al. (2014); Xu et al. (2015). Besides, other types of attention mechanisms are also proposed toward applications in sequence modeling and natural language processing (NLP) such as Cheng et al. (2016); Vaswani et al. (2017).

heading: 4 EXPERIMENTS
text: We evaluate the performance of the imputation-free prediction models trained by GIL against the existing imputation-prediction paradigm on both benchmark and real-world datasets, where the imputation stage employs both state-of-the-art and classic imputation algorithms (IAs), with the details introduced in Section 4.1. We also consider variations of the GIL method proposed in Section 2 to illustrate its robustness and flexibility. The datasets we use include i) MIMIC-III (Johnson et al., 2016) that consists of real-world EHRs obtained from intensive care units (ICUs), ii) a de-identified ophthalmic patient dataset obtained from an eye center in North America, and iii) hand-written digits MNIST (LeCun & Cortes). We also tested on a smaller scaled ICU time-series from 2012 Physionet challenge (Silva et al., 2012) and these results can be found in Appendix D.4. Some of the data are augmented with additional missingness to ensure sufficient missing rates and the datasets we use cover all types of missing data, i.e., missing complete at random (MCAR), MAR and MNAR (Little & Rubin, 2019). We found that the proposed method not only outperforms the baselines on all three datasets under various experimental settings, but also offers better feature representations as shown in Figure 3. We start with an overview of the baseline methods, in the following sub-section, and then proceed to present our experimental findings.

heading: 4.1 VARIANTS OF GIL AND BASELINES
text: Our method (GIL) trains MLPs or LSTMs (depending on the type of data) to directly output the predictions given incomplete inputs without imputation. Following from Section 2.4, we also test on
binary classification tasks using GIL-D which includes the distributional difference term D(ζ+, ζ−), captured by mean squared errors, into the reward function. For baselines, we start with another variant of GIL that uses a simple heuristic – the missing indicator m – to replace the importance a in GIL; we denote it as GIL-H. This helps analyze the advantages of training the models using the importance obtained from the RL policy learned by GIL, versus a heuristic that simply discards the gradients produced by the subset of input entries that are missing.
For other baselines, following the imputation-prediction framework, the imputation stage employs state-of-the-art IAs including MIWAE (Mattei & Frellsen, 2019) which imputes the missing data by training variational auto-encoders (VAEs) to maximize an ELBO, GP-VAE (Fortuin et al., 2020) that couples VAEs with Gaussian processes (GPs) to consider the temporal correlation in time-series, and the generative adversarial network-based method GAIN (Yoon et al., 2018). Further, some classical imputation methods are also considered, including MICE, missForest (MF), k-nearest neighbor (kNN), expectation-maximization (EM), mean-imputation (Mean), zero-imputation (Zero) and carrying-forward (CF). The imputed data from these baselines are fed into the prediction models that share the same architecture as the ones trained by GIL, for fair comparison. For time-series, we also compare to BRITS (Cao et al., 2018) which trains bidirectional LSTMs end-to-end by masking out additional observed data and jointly imposing imputation and classification objectives.

heading: 4.2 MIMIC-III
text: We consider the early prediction of septic shock using the MIMIC-III dataset, following the frameworks provided in recent data-driven works for data pre-processing (Sheetrit et al., 2017; Fleuren et al., 2020; Khoshnevisan et al., 2020). Specifically, we use 14 commonly-utilized vital signs and lab results over a 2-hour observation window to predict whether a patient will develop septic shock in the next 4-hour window. The resulting dataset used for training and testing contains 2,166 sequences each of which has length between 3 and 125, and the overall missing rate is 71.63%; the missing data in this dataset can be considered as a mixture of MCAR, MAR and MNAR. More details can be found in Appendix D.1.
To evaluate performance, we consider the prediction model constituted by a 1024-unit LSTM layer for encoding followed by a dense output layer for inference. Considering that some of the baseline methods are not designed to capture the temporal correlations within the sequences which do not follow a fixed horizon, besides testing with the varied-length sequences (Var-l.) obtained in above, we also test with a fixed-length version (Fix-l.) where the maximum time step for each sequence is set to be 8 (see Appendix D.1 for details). The accuracy1 and AUC for the two tests are summarized in Table 1. When varied-length sequences are considered, GIL(-D) slightly outperforms GAIN and MICE while it significantly dominates the other baseline methods. However, when fixedlength sequences are considered, GAIN and MICE’s performance decreases dramatically and are significantly dominated by that of GIL(-D). This could be caused by omitting the records beyond the 8-th time step as both models are flexible enough to capture such information from that point onwards, while in contrast, the performance increases for kNN, MF and EM which are in general based on much simpler models. On the other hand, the models trained by GIL(-D) was not significantly affected by the information loss. In fact, it emphasizes that applying importance to the gradients during training can enable the model to capture the information behind missing data. Moreover, the dramatically increased performance from GIL-H to GIL(-D) in both tests underscores the significance of training the models using the importance determined by the RL policy learned by GIL(-D), instead of a pre-defined heuristic. Note that according to a recent survey of septic shock predictions (Fleuren et al., 2020), the overall maximum AUC attained by existing data-driven methods, which uses domain knowledge to design models specifically for ICU time-series, is 0.96 and it is comparable to that attained by our method. Finally, GP-VAE and BRITS require the input sequences to have the same horizon, so we only report for the Fix-l. case where they under-perform. These are possibly due to
1All accuracy values in this work are obtained using a decision threshold of 0.5.
the high variation in GPVAE’s uncertainty estimation component, and the mechanism of imposing additional missingness in BRITS2.

heading: 4.3 OPHTHALMIC DATA
text: We consider identifying diabetic retinopathy (DR), an eye disease caused by diabetes, using the de-identified data obtained from Duke Eye Center, constituted by a mixture of image feature vectors and tabular information. Specifically, the data collected from each subject is formulated as a vector with 4,101 dimensions containing two 2,048-length feature vectors from two retinal imaging modalities, optic coherence tomography (OCT) and color fundus photography (CFP), followed by a 5-dimensional data characterizing demographic (e.g., age) and diabetic information (e.g., years of diabetes). Additional details of the dataset and examples of the OCT and CFP images are provided in the Appendix D.2. At most one of the two types of retinal images could be missing due to MAR as sometimes the diagnosis can be performed with a single modality (Cui et al., 2021), while the demographic/diabetic information can be missing due to MCAR. A total of 1,148 subjects are included in this dataset with a 17% missing rate (M.R.)3. We apply additional random masks over both the input image features (by assuming one or two entire image feature vectors are missing) and demographic/diabetic information (by assuming some entries are missing) to increase the missing rate to 25% and 35%, where for each case 3 different random masks are applied.
For this experiment we consider an MLP with 2 hidden layers with 1,000 nodes each for prediction. The results are shown in Table 2 where the mean and standard deviations (in subscripts) of accuracy and AUC are reported. Although GIL and GP-VAE result in similar performance in the 25% missing rate case, GIL still achieves higher accuracy and AUC; these are significantly improved, with significantly lower standard deviation over GP-VAE, when the missing rate increases to 35%. In general, GIL-D outperforms all the other methods, with an exception of the higher standard deviation of accuracy in the 35% M.R. case, which benefits from the flexibility of the framework we proposed. With increased missingness, the performance of most baselines drops significantly and are close to the zero-imputation baseline as the image feature vectors occupy most dimensions of the inputs while imputing with zeros would be a reasonable heuristic.

heading: 4.4 MNIST
text: This work focuses mostly on tabular inputs or time-series; however, we test on MNIST images given its simple structure, which could be classified using MLPs and more importantly, because results from these data are easier to interpret. Specifically, we test on the MCAR version of MNIST where a pre-determined portion of pixels (i.e., out of 50%, 70% and 90%) are masked off uniformly at random from the original images, and the MAR version, which assumes part of the image is always observable and the positions of missing pixels are determined by a distribution conditioned on features extracted from the observable portion following Mattei & Frellsen (2019). For each test, the masks are applied
2The hyper-parameters used to train their models can be found in Appendix D.1 3When calculating missing rates, each feature vector is only counted as a single entry regardless of size.
with 5 different random seeds, resulting in the standard deviations reported in Table 3 which also summarizes obtained results.
We consider using MLPs constituted by two hidden dense layers with 500 nodes each as the prediction models4. Note that in the MCAR setting, GP-VAE achieves slightly higher average accuracy than our method when 50% and 70% of the pixels missing, because it relies on convolutional encoding layers, which intrinsically have advantages on this dataset. However, GP-VAE’s performance is associated with higher standard deviation in both cases and is significantly outperformed by our method when 90% pixels are missing. The models trained by GIL also significantly outperforms the baselines following the imputation-prediction framework. This suggests that most of the IAs fail to reconstruct the input data when the missing rate is high (see Figure 1). Moreover, the errors produced during the imputation are then propagated into the prediction step which results in the inferior performance. In the MAR setting, the task becomes more challenging as entire rows of pixels could be missing. Finally, the model trained by GIL outperforms most of the baselines depending on IAs with an exception of MIWAE, which only slightly outperforms it. Note that GIL-H’s performance is close to GIL in most of the settings due to the simple structure of the MNIST digits where the missing indicator m could be a good estimation of the importance a.
Feature Space Visualization In Figure 3, we use t-distributed stochastic neighbor embedding (tSNE) to visualize the features learned by the prediction model trained by GIL compared to MIWAE and GAIN. We observe that GIL results in more expressive features than the others as it generates clusters with clearer boundaries.

heading: 4.5 CORRELATION BETWEEN IMPUTATION AND PREDICTION PERFORMANCE
text: In this section we study if imputation error is correlated with the downstream prediction performance, under the imputation-prediction framework. Each column of Table 4 is obtained by calculating the Pearson’s correlation coefficient (Freedman et al., 2007) and its p-value, between the imputation mean squared error (MSE) and prediction accuracy, across different imputation methods under the same dataset. It can be observed that the imputation MSE is negatively correlated with prediction performance for most of the datasets and data with higher M.R. tends to have higher degree of negative correlation, which indicates that imputation error could be propagated to downstream tasks.

heading: 5 CONCLUSION
text: We have developed the GIL method, training DL models (MLPs and LSTMs) to directly perform inference with incomplete data, without the need for imputation. Existing methods addressing the problem of missing data mostly follow the two-step (imputation-prediction) framework. However, the error produced during imputation can be propagated into subsequent tasks, especially when the data have high missingness rates or small sample sizes. GIL circumvents these issues by applying importance weighting to the gradients to leverage the information underlying the missingness patterns in the data. We have evaluated our method by comparing it to the state-of-the-art baselines using IAs on two real-world datasets and one benchmark dataset.
4The performance of MIWAE reported here is different from Mattei & Frellsen (2019) because they used a CNN classifier.

heading: ACKNOWLEDGMENTS
text: This research was supported in part by the ONR under agreements N00014-17-1-2504 and N00014-201-2745, AFOSR under award number FA9550-19-1-0169, NIH under NIH/NIDDK R01-DK123062 and NIH/NINDS 1R61NS120246 awards, DARPA, as well as the NSF under CNS-1652544 award and the National AI Institute for Edge Computing Leveraging Next Generation Wireless Networks, Grant CNS-2112562.
We would also like to thank Ruiyi Zhang (Adobe Research) for advice that leads to improved paper presentations, and Ge Gao (North Carolina State University) for sharing insights toward processing the MIMIC-III dataset as well as the use of some imputation baselines.

heading: A EXTENDING GIL TO LSTMS
text: LSTM Model. We also consider the use of an LSTM as the encoder for sequential inputs with varied lengths. Specifically, in this case, we can define Xj = (x>j,1,x > j,2, . . . ,x > j,Tj
)> ∈ RTj×d along with Mj = (m>j,1,m > j,2, . . . ,m > j,Ti
) ∈ {0, 1}Tj×d such that each xj,t ∈ Rd and mj,t ∈ {0, 1}d, where i ∈ [1, N ] and t ∈ [1, Tj ]. Recall that the forward pass of an LSTM follows
ot = σ(Woxj,t + Uoht−1), it = σ(Wixj,t + Uiht−1), ft = σ(Wfxj,t + Ufht−1), gt = tanh(Wgxj,t + Ught−1), ct = ft ct−1 + it gt, ht = ot tanh(ct), (6)
where xt ∈ Rd is the input at time step t, {Wo,Wi,Wg,Wf , Uo,Ui,Ug,Uf} are the weights, σ(x) = 1/(1 + e−x) is the sigmoid function, and is the element-wise product. The output layer is appended after ht to constitute the inference layer, i.e.,
ŷt = φout(Woutht), (7)
where φout is the activation function, Winf = Wout are the weights and t ∈ [1, Tj ]. The following proposition shows that given a loss function E(ŷ,y), the gradients of E w.r.t the encoding weights Wenc = {Wo,Wi,Wg,Wf} can be written in the form of outer products. Note that the gradients for (autoregressive) parameters depending on previous hidden states ht−1, i.e.,, {Uo,Ui,Ug,Uf}, like the bias terms, cannot be written as outer products; thus, these weights are updated following the regular SGD, without importance weighting. The proof is provided following the proposition.
Proposition 2 Given an LSTM as described in (6)-(7) and a smooth loss function E(ŷ,y), the gradients of E w.r.t. Wenc = {Wo,Wi,Wg,Wf} at time t can be written in outer product forms – i.e., it holds that ∂E∂Wo ∣∣∣ t = ∆ot x > t , ∂E ∂Wi ∣∣∣ t = ∆itx > t , ∂E ∂Wg ∣∣∣ t = ∆gt x > t , ∂E ∂Wf ∣∣∣ t = ∆ftx > t , where ∆ot ,∆ i t,∆ g t ,∆ f t are the gradients propagated from the inference layers defined in below, and xt represents the observation obtained at time step t from any sequence Xj ∈ X in general.
Then it follows that if we define x as the t-th observation from any sequence Xj ∈ RTj×d, regardless of it index, (4) can be directly used for training the LSTM encoding weights Wenc = {Wo,Wi,Wg,Wf}. Now we prove the Proposition above.
We first define
∆ot =∂E/∂ht tanh(ct) ot (1− ot), (8) ∆it =∂E/∂ht ot (1− tanh
2(ct)) gt it (1− it), (9) ∆gt =∂E/∂ht ot (1− tanh
2(ct)) it (1− g2t ), (10) ∆ft =∂E/∂ht ot (1− tanh 2(ct)) ct−1 ft (1− ft). (11)
Now to prove the proposition we start from the derivative of the smooth loss function E w.r.t. x, which can be derived as
∂E ∂ht =
[ ∂(Woutht)
∂ht ]> · (12)[
∂E ∂ŷt φ′(Woutht)
] (13)
=W>out
[ ∂E
∂ŷt φ′(Woutht)
] . (14)
Then the derivatives of E w.r.t. the ot, ct, ft, ct−1, it and gt are
∂E ∂ot = ∂E ∂ht ∂ht ∂ot
= ∂E
∂ht tanh(ct) (15)
∂E ∂ct = ∂E ∂ht ∂ht ∂ct
= ∂E
∂ht ot (1− tanh2(ct)) (16)
∂E ∂ft = ∂E ∂ht ∂ht ∂ct ∂ct ∂ft
= ∂E
∂ht ot (1− tanh2(ct)) ct−1 (17)
∂E ∂ct−1 = ∂E ∂ht ∂ht ∂ct ∂ct ∂ct−1
= ∂E
∂ht ot (1− tanh2(ct)) ft (18)
∂E ∂it = ∂E ∂ht ∂ht ∂ct ∂ct ∂it
= ∂E
∂ht ot (1− tanh2(ct)) gt (19)
∂E ∂gt = ∂E ∂ht ∂ht ∂ct ∂ct ∂gt
= ∂E
∂ht ot (1− tanh2(ct)) it. (20)
Now, the derivatives of E w.r.t the weights Wo, Wi, Wg and Wf are
∂E
∂Wo
∣∣∣ t = ∂E
∂ht ∂ht ∂ot ∂ot ∂Wo
(21)
=
[ ∂E
∂ht tanh(ct) ot (1− ot)
] x>t (22)
=∆ot x > t (23)
∂E
∂Wi
∣∣∣ t = ∂E
∂ht ∂ht ∂ct ∂ct ∂it ∂it ∂Wi
(24)
= [ ∂E ∂ht ot (1− tanh2(ct)) gt it (1− it) ] x>t (25) =∆itx > t (26)
∂E
∂Wg
∣∣∣ t = ∂E
∂ht ∂ht ∂ct ∂ct ∂gt ∂gt ∂Wg
(27)
= [ ∂E ∂ht ot (1− tanh2(ct)) it (1− g2t ) ] x>t (28) =∆gt x > t (29)
∂E
∂Wf
∣∣∣ t = ∂E
∂ht ∂ht ∂ct ∂ct ∂ft ∂ft ∂Wf
(30)
= [ ∂E ∂ht ot (1− tanh2(ct)) ct−1 ft (1− ft) ] x>t (31) =∆ftx > t . (32)
Note that since σ′(·) = σ(·)(1− σ(·)), in the transition between (21) and (22) it follows that
∂ot ∂Wo =σ′(Woxt + Uoht−1)x > t (33)
= [ σ(Woxt + Uoht−1) (1− σ(Woxt + Uoht−1)) ] x>t (34)
= [ot (1− ot)] x>t . (35)
Similarly, the transition between (24) and (25) follows
∂it ∂Wi = it (1− it). (36)
At last, the transition between (30) and (31) follows
∂ft ∂Wf = ft (1− ft). (37)
Furthemore, since tanh′(·) = 1− tanh2(·), the transition between (27) and (28) follows
∂gt ∂Wg = tanh′(Wgxt + Ught−1)x > t (38)
=(1− tanh2(Wgxt + Ught−1))x>t (39) =(1− g2t )x>t , (40)
which concludes the proof.

heading: B DESCRIPTION OF ALGORITHM 1
text: The algorithm takes as input the training dataset X , the training targets Y , the missing indicatorsM, the weights of the encoding layers Wenc, the weights for the inference layers Winf , the actor πθ, the critic Qν , learning rates {α, αθ, αν} and training loss function E. Our approach starts by initializing all the parameters Wenc, Winf , πθ, Qν , sampling x ∈ X with the corresponding m ∈ M that will be used for training in the first iteration, obtaining the feature ζ and the prediction ŷ, which constitute the initial state as s = (x,m, ζ, ŷ). In each iteration, first the importance is generated from a behavioral policy β that is conditioned on the target policy πθ, such as the noisy exploration policy proposed in Lillicrap et al. (2015). Then the encoding layer is trained following (4), while the inference layers are trained following the regular gradient descent. After training, the new prediction is obtained following the updated weights and its value is assigned to ŷ, which is then used to generate the reward following the reward function R. Then the training sample x′ for the next iteration is sampled and the corresponding m′, ζ ′, ŷ′ are obtained to constitute the next state s′. Finally, the actor πθ and critic Qν are updated following (5). We refer to Silver et al. (2014); Lillicrap et al. (2015) for more details on actor-critic RL.
C IMPORTANCE VS ATTENTIONS
We now illustrate the connections and distinctions between the importance in the GIL and visual attentions Mnih et al. (2014); Ba et al. (2014); Xu et al. (2015), which are commonly used to train CNNs to focus on specific dimensions of the inputs that are most helpful for making predictions in the context of image captioning and multi-object recognition. In visual attentions, an input image is first encoded into a set of vectors I = {l1, . . . , lL} where each li is a feature vector corresponding to a specific region in the image as they are retrieved from lower convolutional layers. Then, the attention αt,i for time step t and region i ∈ [1, L] is generated following Bahdanau et al. (2014) as αt,i = exp(et,i)/ ∑L j=1 exp (et,k), where et,i = fatt(li,ht−1), fatt is usually an MLP (or the socalled attention model) and ht−1 is the hidden state of a recurrent neural network (RNN) Hochreiter & Schmidhuber (1997) used for prediction. Then, a weighted average of the feature vectors ∑ i lt,ili is fed into the prediction network for further inference where lt,i is a random variable parametrized by αt,i following p(lt,i = 1|lj<t, I) = αt,i, and lj<t represents the historical values of lt,i for all i ∈ [1, L]. Given the definition of the multinoulli variable lt,i, the model is not smooth and thus cannot be trained following the regular back-propagation. So the prediction network is trained to maximize the evidence lower bound (ELBO) which updates the prediction network parameter θ using
∂J ∂θ ≈ 1 n n∑ i=1 [∂ log p(y|l̃i, I) ∂θ + (log p(y|l̃i, I)− b) ∂ log p(l̃i|I) ∂θ ] , (41)
where l̃i are the Monte-Carlo samples drawn from p(lt,i|lj<t, I) and b represents the performance from a baseline E[log p(y|θbaseline)]. It is notable that (41) is equivalent to the REINFOCE algorithm in RL Williams (1992). Specifically, the search space can be interpreted as an MDP which is usually used to model the environment in RL problems – i.e., the state space is constituted by I, the RNN hidden state ht and the historical visitations lj<t; the actions are multinouli variables lt,i; the reward function is defined as the marginal log-likelihood log p(y|lt,i, I); and the policy is characterized by p(lt,i|lj<t, I) Xu et al. (2015); Mnih et al. (2014); Ba et al. (2014). However, these methods cannot be applied directly to our problem as it requires the features I to exclusively associate with specific parts of the inputs spatially, which is attainable by using convolutional encoders with image inputs but intractable with tabular inputs considered in our case (see Section 2.1). Instead, our approach overcomes this issue by directly applying importance, generated by RL, into the gradient space during back-propagation. Moreover, our method does not require to formulate the learning objective as ELBOs and as a result GIL can adopt any state-of-the-art RL algorithm – not limited to REINFORCE only as in Mnih et al. (2014); Ba et al. (2014); Xu et al. (2015).

heading: D EXPERIMENTAL DETAILS AND ADDITIONAL EXPERIMENTS
text: The case studies are run on a work station with three Nvidia Quadro RTX 6000 GPUs with 24GB of memory for each. We use Tensorflow to implement the models and training algorithms. To train the imputation-free prediction models using GIL, we perform a grid search for the model learning rate α ∈ {0.001, 0.0007, 0.0005, 0.0003, 0.0001, 0.00005, 0.00001}, the exponential decay step for α is selected from {1000, 750, 500} and the exponential decay rate for α is selected from {0.95, 0.9, 0.85, 0.8}. The actor πθ and critic Qν in the GIL (i.e., Alg. 1) are trained using deep deterministic policy gradient (DDPG) Lillicrap et al. (2015) where the discounting factor γ = 0.99. We choose the behavioural policy β in GIL as
β(s) = { πθ(s) with probability p1, missing indicator m with probability p2, random action with probability p3.
(42)
The learning rates of the actor απ and the critic αQ are selected by performing grid search from {0.0005, 0.0001, 0.00005, 0.00001} and {0.001, 0.0005, 0.0001} respectively. To train the imputation-free models using GIL-H/GIL-D, we follow the same grid search for α along with its decay steps and decay rates. From the implementation perspective, we replace the missing entries, in the inputs x, with a placeholder value before feeding them into the model. This can avoid value errors thrown by Tensorflow if the input vectors contain NaN’s. However, note these values will not be used to update the parameters. For the state-of-the-art baselines – GAIN5, MIWAE6, GPVAE7 and BRITS8 – we use the implementations published on the Github by the authors. Adam optimizer is used to train all the prediction models for baselines, or to compute (∂E/∂W)SGD for GIL, GIL-H and GIL-D. All the models are trained using a batch size of 128. The details of selecting the other hyper-parameters of each case study can be found in the corresponding sub-section below.
D.1 MIMIC-III
Dataset Formulation The MIMIC-III9 contains EHRs obtained from roughly 40,000 patients who stayed in the ICUs in the Beth Israel Deaconess Medical Center between 2001 and 2012 (Johnson et al., 2016) and septic shock is a severe type of sepsis that results in above 40% mortality rate. Figure 4 illustrates the 2-hour observation window, 4-hour early prediction window (EPW) and the relation between inputs and targets for predictions. Note that each work related to septic shock predictions adopts a slightly different strategy in terms of selecting lab/vital attributes and the lengths of the observation and EPW (Mao et al., 2018; Darwiche & Mukherjee, 2018; Yee et al., 2019; Sheetrit et al., 2017; Liu et al., 2019; Khoshnevisan et al., 2020; Fleuren et al., 2020). In this work, the 4-hour EPW allows including sufficient number of subjects which can ensure statistical significance of the results, while the smaller observation window keeps the task challenging. To formulate the training and testing dataset, we first selected 14 commonly used attributes for predictions as suggested in previous works (Sheetrit et al., 2017; Fleuren et al., 2020; Khoshnevisan et al., 2020), which are constituted by vital signs including temperature, respiratory rate, heart rate, systolic blood pressure,
5https://github.com/jsyoon0823/GAIN 6https://github.com/pamattei/miwae 7https://github.com/ratschlab/GP-VAE 8https://github.com/caow13/BRITS 9Data acquired from https://mimic.mit.edu. Access to this dataset follows the MIT-CLP license.
mean arterial pressure, peripheral oxygen saturation (or SpO2), fraction of inspired oxygen (or FIO2), and lab tests including white blood cell count, serum lactate level, platelet count, creatinine, bilirubin, bandemia of white blood cells, blood urea nitrogen. To form the septic shock cohort, we first identify the patients who were diagnosed with sepsis by indexing with the International Classification of Diseases 9 (ICD-9) code 995.91 and 995.92. Then the patients with septic shock history are selected using ICD-9 code 785.52 with their septic shock onset time determined following the third international consensus for sepsis and septic shock (Sepsis-3) standard (Singer et al., 2016). Finally, we remove the patients whose septic shock onset time was less than 6 hours after admission to the ICUs since the data is not enough to be formulated as sequences pertaining to the 2-hour observation and 4-hour prediction window. Consequently, a total of 1,083 septic shock patients are identified. To form the non-shock (or control) cohort, first 1,083 patients are randomly sampled from all admissions who have at least 2 hours of records excluding the ones who are in the septic shock cohort, then a 2-hour time frame is randomly selected to be used as the observation window. All patients selected following the above procedure are split into 8:2 to formulate the training and testing datasets, which are referred to as the varied-length (Var-l.) sequences in Section 4.2. In the test set, the number of subjects associated with positive and negative labels, respectively, are selected to be equivalent.
Formulation of Fix-l. Sequences To formulate the fixed-length (Fix-l.) sequences, we follow this protocol – if sequence length is less than 8, we pad it with the latest recording available until the length reaches 8, otherwise we truncate the length to 8 by discarding the data from that point on-wards. We specifically choose this threshold because most of the sequences have length ≤ 8, as shown in Figure 5.
Type of Missingness The missing data in this dataset can be considered as a mixture of MCAR, MAR and MNAR. Specifically, recordings from human mistakes or malfunctioned equipment manifest as MCAR, the dependency across vital signs and lab results give rise to MAR (e.g., specific lab tests are ordered only if abnormalities found in the related vital readings or other lab results), and the mismatch among the sampling frequency of some periodically recorded attributes such as temperature (obtained hourly) and blood test (conducted daily) can be considered as MNAR (Little & Rubin, 2019; Mamandipoor et al., 2019).
Training Details To train the prediction models for both GIL and baselines, we use maximum training steps of 2,000 and 4,000 for varied-length and fixed-length sequences respectively. For imputation of the missing values, we train the imputation method MIWAE using the learning rates {0.001, 0.0001} with other hyper-parameters provided by the authors in the code. To train GAIN, we use the default learning rates provided by the authors to train the generator and discriminator. GAIN contains one hyper-parameter that balances between the two losses LG and LM defined in Yoon et al. (2018) which is selected from {0.1, 1, 10, 100}. To train GP-VAE on the Fix-l. case, the set of hyperparameters we use contain latent_dim = {6, 12, 35}, encoder_size = 256, 256, decoder_size = 256, 256, 256, window_size = {3, 4, 6}, beta = {0.2, 0.5, 0.8}, sigma = 1.005, length_scale =
{2, 4, 7}. To train BRITS we used the recommended impute_weight = 0.3, label_weight = 1.0. The term D(ζ+, ζ−) included in the reward function for GIL-D is defined as in the follows. Assume that in each training epoch b/2 inputs with label 0 and b/2 inputs with label 1 are sampled from the dataset, where b is the batch size. We use F+ = (ζ+1 , . . . , ζ + b/2) ∈ R
b/2×e to denote the features associated with positive labels (i.e., 1) in the current batch and F− = (ζ−1 , . . . , ζ − b/2) ∈ R
b/2×e are , where e is the dimension of each individual feature ζ. Then we define D as
MSE(F+[0 : b/4], F−[0 : b/4]) +MSE(F+[b/4 : b/2], F−[b/4 : b/2])
−MSE(F+[0 : b/4], F+[b/4 : b/2])−MSE(F−[0 : b/4], F−[b/4 : b/2]), (43)
where the slicing index follows the syntax from Python, e.g., F+[0 : b/4] corresponds to the 0-th to (b/4− 1)-th rows of F+.
D.2 OPHTHALMIC DATA
Dataset Introduction This dataset is originally constituted by OCT images, CFP images and patient EHRs including demographic information and test results related to diabetes. A total of 1,148 subjects are included in this dataset. OCT refers to the two- or three-dimensional images capturing the retinal architectures of the eyes that are scanned by low-coherence lights. The 3D OCT image is usually called the volume scan and it is constituted by 61 or 101 2D image slices depending on the spec of the scanner. Examples of 3D and 2D OCT scans are shown in Figure 6. In this experiment the volume scans we use contain the 61 slices, while we only consider the center slice (31th) 2D image. The CFP images are the color images captured by a fundus camera showing the condition of the interior surface of the eye, where an example is shown in Figure 7.
Dataset Formulation To formulate the data into tabular inputs, the 2D OCT and CFP retinal images are first fed in to two inception-v3 (Szegedy et al., 2016) CNNs, pre-trained with similar images provided by (Kermany et al., 2018) and (Kaggle, 2015) respectively, and the image feature vectors with dimension of 2,048 each are output by the global average pooling layer. Then the patient EHR information include age (integer), sex (boolean), length of diabetic history (integer), A1C result (integer), which measures blood sugar levels, and if insulin has been used (boolean) constitute a 5-dimensional vector that is concatenated to the end of the OCT and CFP feature vectors. We split all the subjects into a training cohort and a testing cohort following a ratio of 9:1. In the test set, the number of subjects associated with positive and negative labels, respectively, are selected to be equivalent. The raw images from the training cohort are augmented through cropping and rotation before feeding into inception-v3. All the data from the testing cohort are not augmented. Figure 7 illustrates how the training and testing inputs are formulated, as well as the input-output relation of the prediction models.
Training Details To train the prediction models for both GIL and baselines, we use maximum training steps of 2,000. For imputation of the missing values, MIWAE is trained using the learning rates {0.001, 0.0001} with other hyper-parameters provided by the authors in the
code. GAIN is trained using hyper-parameters selected from {0.1, 1, 10, 100}. GP-VAE is trained by applying the arguments latent_dim = 256, encoder_size = 256, 256, decoder_size = 256, 256, 256, window_size = 3, beta = 0.8, sigma = 1 to its code published in Github (https://github.com/ratschlab/GP-VAE). MF, EM and kNN are not included in this experiment due to the high dimensionality of the inputs which makes these algorithms very computational inefficient as imputation results cannot be produced within days.
D.3 MNIST
Training Details To train the prediction models for both GIL and baselines, we use maximum training steps of 1,0000. For imputation of the missing values, MIWAE is trained using the learning rates {0.001, 0.0001} with other hyper-parameters provided by the authors in the code. GAIN is trained using hyper-parameters selected from {0.1, 1, 10, 100}. GP-VAE is trained by applying the arguments latent_dim = 256, encoder_size = 256, 256, decoder_size = 256, 256, 256, window_size = 3, beta = 0.8, sigma = 1 to its code published in Github (https://github.com/ratschlab/GP-VAE). MF, EM and kNN are not included in this experiment due to the high dimensionality of the inputs which makes these algorithms very computational inefficient as imputation results cannot be produced within days.
Additional Results We also trained GIL-D on the MCAR version of MNIST dataset and the performance are shown in Table 5 below.
Other than MIMIC-III, we also tested on a smaller scaled ICU time-series from 2012 Physionet challenge Silva et al. (2012) which contain data obtained from 12,000 patients. We use the data pre-processed and open-sourced by Fortuin et al. (2020).10 For each patient the values of 35 different attributes (e.g., blood pressure, temperature) are recorded over a 48-hour window. As a result, the data for each patient can be formulated into a matrix in R48×35, i.e., the sequence length across patients are the same. This dataset has an overall 78.5% missing rate and a binary label is assigned to each patient where 87% of them are 0’s and 13% are 1’s. Therefore we include average precision (AP) which is calculated from the precision-recall curve as an additional metric to evaluate the performance toward imbalanced labels. For this dataset, we consider a 1024-unit LSTM layer for encoding and a dense output layer for inference. It can be observed from Table 6 that although our method achieves the highest accuracy and 2-nd highest AP, while most methods perform very close to mean- and zeroimputation. This could be caused by the its simple structure, as all the sequences share the same time horizon, which significantly reduces the difficulties for the classification task. Moreover, the highly imbalanced labels and the small population result in the performance to be hardly distinguishable across different methods. And this is the reason that we focus on the MIMIC-III dataset to study the strengths and shortcomings of the methods.
10https://github.com/ratschlab/GP-VAE

heading: E PROOF FOR PROPOSITION 1
text: Here we prove a generalized version of Proposition 1 as stated in the following proposition, which shows that the gradient of the loss function E w.r.t. any layer in the MLP can be written in the outer product format.
Proposition 3 Given an MLP as in (1) and a smooth loss function E(ŷ|y), the gradients for any hidden layer Wi, i ∈ [1, k] can be represented as
∂E
∂Wi = ∆iq
> i−1, (44)
where qi−1 = φi−1(Wi−1qi−2) is the output from the i − 1-th layer with q0 = x, and ∆i = W>i+1∆i+1 φ′i(Wiqi−1) with ∆out = ∂E∂ŷ φ ′ out(Woutq2k) and denotes the element-wise (Hadamard) product.
Now we prove Proposition 3.
We consider the MLP characterized by
ŷ = φout(Woutφk(Wk . . . φ2(W2φ1(W1x)))) (45)
where x ∈ Rd represents the input to the model, Wi is the weight matrix for the i-th hidden layer, and φi is the activation function of the i-th hidden layer.
We start with deriving the derivatives for the output layer
∂E ∂Wout = ∂E ∂ŷ ∂ŷ ∂Wout (46)
=[ ∂E
∂ŷ φ′out(Woutqk)] · q>k (47)
=∆out · q>k , (48)
where qk is the output from the k-th hidden layer.
Now we show the derivative of E w.r.t. the k-th hidden layer
∂E ∂Wk = ∂E ∂ŷ ∂ŷ ∂qk ∂qk ∂Wk
(49)
=W>out[ ∂E
∂ŷ φ′out(Woutqk)] · ∂qk ∂Wk
(50)
=W>out∆out · ∂qk ∂Wk
(51)
=[(W>out∆out) φ′k(Wkqk−1)] · q>k−1 (52) =∆k · q>k−1 (53)
We now prove Proposition 3 by induction. First assume that for j ∈ [2, k], the derivative of E w.r.t. the j-th hidden layer is
∂E
∂Wj =[(W>j+1∆j+1) φ′j(Wjqj−1)] · q>j−1 (54)
=∆j · q>j−1. (55)
Moreover, let us assume that
∂E ∂qj =W>j+1[ ∂E ∂qj+1 φ′j+1(Wj+1qj)] (56)
=W>j+1∆j+1. (57)
Now we need to show that for j − 1, it holds that ∂E
∂Wj−1 =[(W>j ∆j) φ′j−1(Wj−1qj−2)] · q>j−2, (58)
∂E
∂qj−1 =W>j [
∂E ∂qj φ′j(Wjqj−1)] (59)
=W>j ∆j . (60)
We first prove that (60) holds, i.e.,
∂E ∂qj−1 = ∂E ∂qj ∂qj ∂qj−1
(61)
=W>j+1∆j+1 ∂qj ∂qj−1
(62)
=W>j [W > j+1∆j+1 φ′j(Wjqj−1)] (63) =W>j ∆j , (64)
which is equivalent to (60).
To prove that (58) holds, we start with
∂E ∂Wj−1 = ∂E ∂qj ∂qj ∂qj−1 ∂qj−1 ∂Wj−1
(65)
=W>j+1∆j+1 ∂qj ∂qj−1 ∂qj−1 ∂Wj−1
(66)
=W>j [W > j+1∆j+1 φ′j(Wjqj−1)] ∂qj−1 ∂Wj−1
(67)
=W>j ∆j ∂qj−1 ∂Wj−1
(68)
=[(W>j ∆j) φ′j−1(Wj−1qj−2)] · q>j−2, (69)
which is equivalent to (58).

heading: F ADDITIONAL RELATED WORKS
text: RL for Optimization. This work is also related to approaches that use RL to solve optimization problems, or to improve existing solvers. For example, Boyan & Moore (2000) proposes a framework to improve the efficiency of local search algorithms by predicting the search outcomes and biasing toward directions that provide higher returns, which is approached by estimating value functions of the MDP corresponding to the optimization problem. Recently, Li & Malik (2016) trains RL agents to perform gradient descent steps over linear regression and neural network models. Similarly, Wei et al. (2020) introduces a policy network to automatically determine parameters for the plug-and-play frameworks (Sreehari et al., 2017) which can solve inverse imaging problems. Note that these methods require complete inputs and cannot be adapted to the problem we consider.

heading: G ADDITIONAL ABLATION STUDY
text: In this section, we compare GIL against an ablation baseline that applies the importance directly to the inputs x, instead of the gradient space. Specifically, suppose a differential parametric function hψ(x), parameterized by ψ, is multiplied with x element-wise. Then, as opposed to (4), the gradient updates w.r.t. the encoding layers can formulated as
W′enc ←Wenc − α∆ · (x> hψ(x)>). (70)
Consequently, all model elements (i.e., hψ, Wenc and Winf ) could be trained by gradient descent, without introducing RL policies.
In Table 7 it shows the performance of the model described above over the ophthalmic dataset, as well as MNIST digits with 70% and 90% missing rate (M.R.). Specifically, hψ(x) is captured by a neural network with the same architecture as the policy πθ in GIL. The training and testing are conducted following the same convention (e.g., random seeds that determine the missing entries and hyper-parameter search) as introduced in Sections 4.3, 4.4 and Appendix D.2, D.3. It can be observed that the ablation baseline attained similar performance as to GIL-H and zero imputation. The reason that leads to these results would be that the gradients for training hψ(x) still follow the outer product format ∆ · x> which are left to be accounted for. Specifically, our method is built on top of the idea, and the experiments in Section 4 also support that if part of the inputs are missing they may not provide sufficient information to train the prediction models, given the outer product format of the gradients. Similarly, the gradients for hψ(x) also follow this format and hψ(x) may not be properly trained directly on incomplete inputs x, which could in general limit the overall performance of the prediction model.
Additional rationales for using RL to obtain gradient importance in GIL. The ablation baseline above is closely related to visual attention (VA) models as discussed in Appendix C, and could be seen as a preliminary version of them where hψ(x) is used to re-weight elements in the inputs x. VA techniques benefit from the fact that features learned by CNNs are spatially correlated to the inputs, so the attentions could be directly applied to learned features. To achieve this, VA formulates the objective as an ELBO which is shown equivalent to the objective of a basic policy gradient RL algorithm, REINFORCE (Mnih et al., 2014). However, it was not clear how such VA methods could be adapted to the problem we consider, as the features learned by MLPs or LSTMs are not spatially correlated with inputs.
More on importance versus attentions. The other baseline, BRITS (Cao et al., 2018), uses bidirectional RNN to process time-series, with attentions applied to the hidden states of RNNs, followed by optimizing over imputation and prediction objectives jointly during training. BRITS is considered as the state-of-the-art method that could predict with incomplete time-series inputs end-to-end. However, we showed that BRITS is outperformed by our method on the MIMIC dataset; see Table 1. The reason would be that BRITS requires masking off part of the observed inputs to constitute the imputation objective; thus, the information provided for model training is even more limited given the intrinsic > 70% M.R. of MIMIC.

Review summaries:
The paper proposes gradient importance learning for predicting labels on incomplete data. The problem is: during both training and testing, certain dimensions of the samples are masked randomly and the masks (i.e. the "missing indicators") are known for the algorithm. Among existing methods in this problem, the paper claims that it is the first "imputation-free" way to solve this problem. The paper focuses on MLP and LSTM and adds an (element-wise) multiplicative parameter matrix on the gradient of the first layer (the one that is mostly closed to the data). The matrix is trained using reinforcement learning. The authors also connect its proposal to visual attention. Experiments are conducted on tabular, time series, and a simple image dataset. 

The paper proposes a method to train on data with missing features. It does so with a single model that can handle missing features, not with extra imputation methods. This is done by weighting the gradient update of the first weight matrix of the neural network with a vector a that is produced by an RL agent. This agent's reward in turn is the performance of the model after the update.

This paper proposes a method to do MLP / LSTM-based inferences/predictions based on incomplete data. While typical methods would impute and then predict, the proposed method predicts directly based on incomplete data. The method learns an “importance” matrix which is multiplied with the first weight matrix of the neural network. The learning is done using reinforcement learning, with the negative prediction error being the reward. This method is compared with new and classical prediction methods for incomplete data on a range of datasets and is shown to be the best most of the time.

This paper targeted at the missing data issue in time series data and proposed a imputation-free method to handle missing data. More specifically, the authors proposed a gradient importance learning method named GIL to weigh the gradients for different parameters using reinforcement learning. Experiments on one tabular dataset and two image datasets demonstrated the effectiveness of the proposed method.