HOW TO DEAL WITH MISSING DATA IN SUPERVISED DEEP LEARNING?
Year: 2021

Abstract:
The issue of missing data in supervised learning has been largely overlooked, especially in the deep learning community. We investigate strategies to adapt neural architectures for handling missing values. Here, we focus on regression and classification problems where the features are assumed to be missing at random. Of particular interest are schemes that allow reusing as-is a neural discriminative architecture. To address supervised deep learning with missing values, we propose to marginalize over missing values in a joint model of covariates and outcomes. Thereby, we leverage both the flexibility of deep generative models to describe the distribution of the covariates and the power of purely discriminative models to make predictions. More precisely, a deep latent variable model can be learned jointly with the discriminative model, using importance-weighted variational inference, essentially using importance sampling to mimick averaging over multiple imputations. In low-capacity regimes, or when the discriminative model has a strong inductive bias, we find that our hybrid generative/discriminative approach generally outperforms single imputations methods.

text: 1 INTRODUCTION
Missing data affect data analysis across a wide range of domains and the sources of missing values span an equally wide range. Recently, deep latent variable models (DLVMs, Kingma & Welling, 2014; Rezende et al., 2014) have been applied to missing data problems in an unsupervised setting (e.g. Rezende et al., 2014; Nazabal et al., 2020; Ma et al., 2018; 2019; Ivanov et al., 2019; Mattei & Frellsen, 2018; 2019; Yoon et al., 2018; Li et al., 2019; Ipsen et al., 2021; Ghalebikesabi et al., 2021), while the supervised setting has not seen the same attention. The progress in the unsupervised setting is focused on inference and imputation and can therefore be useful as an imputation step before learning a discriminative model. Traditionally the focus has also been on inference and imputation either as a goal in itself or before supervised learning on the (possibly multiple) imputations and observed data (Rubin, 1976; 1996; Little & Rubin, 2019).
Supervised learning in the presence of missing values has different goals and pose different challenges than inference and imputation (Josse et al., 2019; Le Morvan et al., 2020a; 2021). Here the overall aim is to minimize an expected prediction error. However, optimal single imputation does not necessarily lead to optimal
prediction in terms of minimizing a prediction error (Josse et al., 2019). One challenge is that predictions based on inputs with missing values can be ambiguous, that is, the conditional distribution of the missing values given the observed may be multimodal and the optimal prediction may change with the mode. With single imputation, the conditional distribution over the missing data is discarded, and the optimal single imputation can no longer reflect this ambiguity, see figure 2. Even the optimal single imputation leads to biased parameter estimates compared to the fully observed case (Bertsimas et al., 2021), and may in some cases lie outside the distribution of the data. Another challenge is that the number of missing value patterns grows combinatorially with the number of features p, requiring 2p submodels to fit the Bayes predictor in the linear case (Le Morvan et al., 2020a). In multiple imputation (Rubin, 1996), several imputations are drawn from the posterior predictive distribution of the missing values, reflecting the full conditional distribution of the missing values and thus the uncertainty about what is missing. This allows for uncertainty estimates in downstream tasks such as prediction. This in turn requires fitting as many discriminative models as the number of imputations.

heading: 1.1 CONTRIBUTIONS
text: In order to address supervised deep learning with missing values we develop the supervised missing data importance-weighted autoencoder (supMIWAE) bound, based on the approximate maximum likehood techniques used by Burda et al. (2016) and Mattei & Frellsen (2019). This is a scalable approach to marginalizing over missing values in a joint model of covariates and outcomes
pφ,θ(y,x obs,xmiss) = pφ(y|xobs,xmiss)pθ(xobs,xmiss). (1) The covariate model pθ(xobs,xmiss) = ∫ p(xobs,xmiss|z)p(z) dz is a DLVM in this work, but can be any probabilistic model imposing a joint density over covariates. The outcome model pφ(y|xobs,xmiss) is any neural discriminative architecture that would have been used in the complete-data case, parameterizing a density over the outcomes. The graphical model is shown in figure 1 along with its computational structure.
Once the joint model has been trained using the supMIWAE bound, the conditional distribution of the outcome given the observed parts of the input can be found as
pφ(y|xobs) = ∫ pφ(y|xobs,xmiss)pθ(xmiss|xobs) dxmiss. (2)
This is approximated using importance sampling techniques, averaging predictions over multiple importance samples (akin to multiple imputations, see appendix C for a deeper discussion) from the generative model.
The model can be trained end-to-end or a pretrained generative model can be coupled with a discriminative model. Having a joint model allows for supervised imputations, where available labels can help guide the imputations from the generative model by adjusting the importance weights accordingly.
Joint models over covariates and labels have previously been used to marginalize over missing values in the covariates, using Gaussian mixture models either as the model over the covariates or directly as the joint model over covariates and outcomes (Ghahramani & Jordan, 1995; Ghahramani & Hinton, 1996; Ahmad & Tresp, 1993; Tresp et al., 1994; 1995; Śmieja et al., 2018). Our approach shows how to use more flexible generative models, such as DLVMs, provides an efficient optimization procedure and allows for keeping any deep discriminative architecture unchanged.

heading: 2 BACKGROUND AND NOTATION
text: We define the random variable x = (x1, . . . , xp) ∈ X which takes values in a p-dimensional feature space X = X1 × · · · × Xp. There is a corresponding (possibly vector valued) response variable y ∈ Y . A missing process obscures parts of x resulting in the mask random variable s ∈ {0, 1}p where
sj = { 1 if xj observed, 0 if xj missing.
(3)
We let obs(s) denote the indices of the non-zero entries of s and miss(s) denote the indices of the zero-entries of s, such that xobs(s) is the set of observed elements of x, and xmiss(s) is the set of missing elements of x. For simplicity we will omit the s and write xobs, xmiss respectively, whenever the context is clear.
Following (Yoon et al., 2018; Le Morvan et al., 2020b; Ghalebikesabi et al., 2021) we define the incomplete random variable x̃ = (x̃1, . . . , x̃p) ∈ X̃ which takes values in X̃ = (X1∪{na})×· · ·× (Xp ∪ {na}) where missing values are represented by the symbol na, such that na · xj = na and na · 0 = 0. We typically only have access to x̃, where
x̃ = x s + na (1− s) (4) and is the Hadamard product. Finally, we are given n i.i.d. copies of the random variables x̃ and y which we collect in a dataset D = {x̃i,yi}ni=1, or alternatively D = { xobsi ,yi }n i=1 .
We make the additional assumption that the data are missing at random (MAR, see e.g. Seaman et al., 2013; Little & Rubin, 2019). Specifically, this means that we assume that s and xmiss are independent given xobs. This assumption allows to avoid to explicitly model the missing mechanism. Our approach could be extended beyond MAR by using a generative model fit for this purpose, like the not-MIWAE of Ipsen et al. (2021), or the deep pattern-set mixture of Ghalebikesabi et al. (2021).

heading: 2.1 CHALLENGES WHEN TRAINING DISRIMINATIVE MODELS WITH MISSING DATA
text: Our predictive model will be defined through a mapping fφ : X → H used to parameterize a conditional distribution pφ(y|x) = Ψ(y|fφ(x)). (5) Here (Ψ(·|η))η∈H is a parametric family of distributions over the outcome space Y , such as a Gaussian distribution in regression or a categorical distribution in classification. Under the maximumlikelihood framework (or equivalently the logarithmic scoring rule), an optimal mapping f∗φ within a class F = (fφ)φ∈Φ parameterized by φ ∈ Φ is defined as
f∗φ ∈ arg max fφ∈F
Ep∗(x,y) [ log Ψ(y|fφ(x)) ] , (6)
where p∗(x,y) is the true data generating distribution. When the data are complete, this is typically approximated by maximizing the log likelihood of the parameters φ given the data Dtrain,
`(φ) = ∑
(x,y)∈Dtrain
log pφ(y|x). (7)
When the covariates have missing values the log-likelihood cannot be maximized directly as the likelihood is not defined since pφ(y|x) depends on the full input vector. The simplest approach is instead to learn separate mappings for each missing pattern. This reduces the amount of data available for the training of each mapping and does not scale well with the dimensionality of the input as in general 2p networks need to be trained. An often used approach is instead to impute the missing values using a single imputation (for instance, using Epθ(xmiss|xobs)[xmiss] or an approximation of it), and then map the concatenation of observed and imputed values to approximate the conditional distribution
pφ(y|xobs) ≈ pφ(y|x = (xobs,Epθ(xmiss|xobs)[x miss])) = Ψ(y|fφ(xobs,Epθ(xmiss|xobs)[x miss])). (8)
Instead of Epθ(xmiss|xobs)[xmiss], which is the optimal imputation under the mean squared error, one could use any kind of imputation, for instance using a constant, or the mean.
This general approach is adequately called impute-then-regress by Bertsimas et al. (2021) and Le Morvan et al. (2021). While it will be Bayes-consistent given a powerful enough classifier (Le Morvan et al., 2021), it leads to biased parameter estimates compared to having complete data (Bertsimas et al., 2021). This is illustrated in figure 2 where the optimal (under the mean squared error) single imputations are placed entirely within one of the classes, requiring a biased decision surface to properly reflect the class label probabilities for the records with missing values. From a probabilistic perspective, computing pφ(y|xobs) requires marginalizing over the missing features as in equation (1), that is
pφ(y|xobs) = Epθ(xmiss|xobs)[pφ(y|x obs,xmiss)]. (9)
Notice the difference to single imputation since, in general
Epθ(xmiss|xobs)[pφ(y|x obs,xmiss)] 6= pφ(y|x = (xobs,Epθ(xmiss|xobs)[x miss])), (10)
except in some pathological cases (e.g. when pφ is linear or pθ is a Dirac distribution). Therefore, if the discriminative model is nonlinear and the generative model is complex enough, then using single imputation will give a very different result than marginalising over the missing features. This is exemplified in figure 2, where a classifier trained using the optimal imputation finds a solution that is very different to what would be found by one trained without missing data. Our technique, marginalizing the missing values, allows to recover the same decision surface as a classifier trained without missing data.

heading: 2.2 OUR JOINT MODEL AND ITS SUPMIWAE VARIATIONAL BOUND
text: As mentioned in the introduction, we posit a joint model over covariates and outcome, that will make use of a latent variable z:
pφ,θ(y,x obs,xmiss, z) = pφ(y|xobs,xmiss)pθ(xobs,xmiss|z)pθ(z). (11)
Under the MAR assumption, the relevant substitute of the likelihood is the likelihood of the observed data
pφ,θ(y,x obs) = ∫ pφ(y|xobs,xmiss)pθ(xobs,xmiss|z)pθ(z) dzdxmiss. (12)
The integral in equation (12) is in all but the simplest cases analytically intractable and direct maximum likelihood methods for learning the parameters (θ,φ) cannot be used. In order to learn the parameters of the joint model, we turn to amortized importance weighted variational inference, where a lower bound is maximized instead of the log likelihood itself (Burda et al., 2016; Domke & Sheldon, 2018). It is based on the fact that unbiased estimates of a likelihood can be turned into lower bounds of the log-likelihood. Indeed, if R(xobs,y) is a random variable such that E [ R(xobs,y) ] = pθ,φ(x obs,y), then E [ logR(xobs,y) ] ≤ logE [ R(xobs,y) ] = log pθ,φ(x
obs,y), which means that E [ logR(xobs,y) ] is a lower bound of the log likelihood. A simple way to find a suitable R(xobs,y) is using importance sampling from a variational distribution over z. More specifically, let
RK(x obs,y) =
1
K K∑ k=1 pφ(y|xobs,xmissk )pθ(xobs|zk)pθ(zk) qγ(zk|xobs) , (13)
where qγ(z|xobs) is the variational distribution (learnable proposal) and (xmissk , zk)k∈{1,...,K} are i.i.d. samples from pθ(xmiss|z)qγ(z|xobs). Then we have that E(xmissk ,zk) [ RK(x obs,y) ] = pθ,φ(x obs,y) and E(xmissk ,zk) [ logRK(x obs,y) ]
is a lower bound on the log likelihood. The variational distribution can be parameterized via a neural network with weights γ ∈ Γ, similarly to Mattei & Frellsen (2019), Nazabal et al. (2020), and Ipsen et al. (2021).
Now that we have a way to bound each per-datum log-likelihood, we can define our supMIWAE variational bound LK(θ,φ,γ) that is a lower-bound of the log-likelihood of the observed data
LK(θ,φ,γ) = n∑ i=1 E(xmissk ,zk) [ logRK(x obs i ,yi) ] ≤ n∑ i=1 log pφ,θ(yi,x obs i ). (14)
The lower bound LK(θ,φ,γ) is a Monte Carlo objective that can be maximized instead of the log likelihood itself, using techniques from stochastic optimization (see e.g. Mohamed et al., 2020). As detailed in appendix D, it benefits from the usual theoretical advantages of importance-weighted variational bounds: LK(θ,φ,γ) converges monotonically to the loglikelihood ∑n i=1 log pφ,θ(yi,x obs i ) at speed 1/K. This means that, the larger the K, the closer our objective will be to the true likelihood.
Note that the contribution of any fully observed data point (x,y) will be simply
log pφ(y|x) + E(zk)
[ log ( 1
K K∑ i=1 pθ(x|zk)pθ(zk) qγ(zk|x)
)] , (15)
which is just the sum of the discriminative likelihood and the standard IWAE bound of Burda et al. (2016). An interesting property of our bound is that the discriminative and generative parts can be trained jointly with a single and statistically sound objective. However, if we were to encounter a data set without missing data, equation (15) means that the training of the two parts of the models would be decoupled: the discriminative model would be the same as a standardly trained discriminative model, and the generative model would be the same as a standard IWAE-trained DLVM.

heading: 2.3 PREDICTION
text: Once the model has been jointly trained using the supMIWAE bound, it is possible to use it to do prediction with missing data. Indeed self-normalized importance sampling can be used to approximate the conditional
pφ(y|xobs) ≈ K∑ i=1 wkpφ(y|xobs,xmissk ), (16)
where
wk = rk
r1 + ...+ rK , and rk = pθ(x obs|zk)pθ(zk) qγ(zk|xobs) , (17)
and (xmissk , zk)k∈{1,...,K} are i.i.d. samples from pθ(x miss|z)qγ(z|xobs). Again, this technique mimicks multiple imputation: several imputations are generated using the generative model, are then fed to the classifier, and the predictions of each imputation are averaged.

heading: 3 RELATED WORK
text: Marginalizing over missing covariates Different approaches have previously been taken to marginalize over missing values in joint models over features and labels. In a series of publications (Ahmad & Tresp, 1993; Tresp et al., 1994; 1995), a closed-form solution using Gaussian Basis Function networks was used to approximate this marginalization, while keeping the neural discriminative model fixed. Ghahramani & Jordan (1995) took a different approach, modeling the joint distribution p(y,x) directly as a mixture model, obtaining p(y|xobs) by conditioning on observed quantities in the joint distribution, thus removing the explicit discriminative neural architecture. In the context of kernel methods, Dick et al. (2008) proposed to learn a distribution over missing values by minimising the regularized empirical risk. Śmieja et al. (2018) trained a Gaussian mixture model (GMM) and discriminative model jointly, and in place of any missing values the activation of the corresponding input neuron was set to the average activation over the GMM conditioned on observed values.
Consistency of single imputation A review of approaches to handling missing data in (non-deep) supervised learning was given by Josse et al. (2019). Here it was shown that under some assumptions about the capacity of the learner, constant imputation is asymptotically consistent in the supervised setting. Le Morvan et al. (2020b) investigated the case of a linear predictor on covariates with missing data, showing that in the presence of missing, the optimal predictor may not be linear and how constant imputation of each feature can be optimized with regards to the model loss. The linear case was further investigated in (Le Morvan et al., 2020a), deriving the analytical form of the optimal predictor and proposing NeuMiss networks to approximate this. Bertsimas et al. (2021) analyzed the impute-then-regress situation and noted that while this paradigm is widely used in practice it incurs a bias in the parameter estimates, but mean imputation is still asymptotically consistent. Le Morvan et al. (2021) showed that the impute-then-regress approach is asymptotically Bayes optimal given a learner with large capacity and that this holds for all missing mechanisms, even missing not at random (Rubin, 1976). However, the choice of imputation function affects the difficulty of the subsequent regression task, when not in the asymptotic regime. As a consequence they proposed to learn imputation and regression jointly, where the missing values are handled by a NeuMiss network.
Other approaches for missing covariates Yi et al. (2019) tackled the issue of sparsity, and specifically large variations in sparsity, by introducing sparsity normalization. This addresses the issue of model output covarying with the sparsity level in the input. Ma et al. (2018; 2019) used a permutation invariant setup to avoid imputing missing data in the input of a variational autoencoder. This approach can be readily extended to the supervised setting, using the permutation invariant setup as a modified input layer as we show in appendix A.3. Some tree approaches to regression/classification can use observations with missing values directly (see e.g. Twala et al., 2008; Chen & Guestrin, 2016; Josse et al., 2019; Gómez-Méndez & Joly, 2021). Sparse coding was used by Caiafa et al. (2020) to train a dictionary and a classifier simultaneously.
DLVMs and missing values Deep latent variable models have been applied to missing data problems in the unsupervised setting, focused on inference, single and multiple imputations. Early works assumed that a DLVM had previously been trained on complete data (Rezende et al., 2014; Mattei & Frellsen, 2018), but a variety of techniques to handle incomplete training sets was then proposed, both in M(C)AR (Yoon et al., 2018; Ma et al., 2018; 2019; Ivanov et al., 2019; Mattei & Frellsen, 2019; Li et al., 2019) and MNAR (Ipsen et al., 2021; Ghalebikesabi et al., 2021) cases. In (Li & Marlin, 2020) irregularly sampled time-series are approached as a missing data problem and they propose VAE and GAN based models for handling this. The proposed models can also be used for supervised learning by training a classifier on the latent representations, either by training the model jointly or using a pretrained encoder. Without missing data, the idea of jointly training a VAE and a discriminative model, has been explored by Ji et al. (2020) or Joy et al. (2021).

heading: 4 EXPERIMENTS
text: We now evaluate discriminative models trained using the supMIWAE bound on a range of supervised learning tasks. Throughout the experiments the generative part of the model is pretrained using the MIWAE bound of Mattei & Frellsen (2019). The pretrained DLVM can be used to draw single imputations, referred to as MIWAE single imputation, or it can be used as the generative part in the supMIWAE computational structure. Here, the generative part of the model is kept fixed while updating the discriminative part of the model according to the supMIWAE bound. As the DLVM in MIWAE single imputations and the supMIWAE are similar, any difference in performance between the two is attributed to their different strategies for handling missing values. Other strategies for handling missing data, used for comparison in this work, are briefly described in appendix A. Missing values are introduced in the training, validation and test sets according to the missing mechanism used in the given experiment. Predictive performance is measured on a test-set with missing values. In all experiments, there is a computational overhead compared to single imputation of training the supMIWAE that scales linear with the number of importance samples K, c.f. equation (13).

heading: 4.1 2D CLASSIFICATION
text: A qualitative analysis of the bias in the learnt parameters of the discriminative model, due to single imputation, is set up on simple 2D datasets, see figure 8 of appendix E. The simplicity allows us
to directly inspect the impact of different imputation strategies on the learnt decision surface and compare this to the marginalization procedure used in the supMIWAE.
All datasets consist of a training set with 3k records and validation and test sets with 1k records. An MCAR missing process is introduced in the horizontal coordinate, where each element becomes missing with probability m = 0.5. Further training details are in appendix E. The datasets are such that there will be inherent ambiguity about the class label for some of the records with missing values. For example in the Burger dataset, all records with a missing value will have p(y = 1|xobs) = 2/3. In figure 3 imputations from different imputation models are shown along with the corresponding learnt decision surface. Partially observed examples are shown in the margins of the plots as sticks, colored by their class label. All methods use the same discriminative architecture, only the imputation strategy differs. For the supMIWAE, multiple imputations are shown for three partially observed records, alongside their kernel density. These multiple imputations reflect the appearance of importance samples used during training, for these three partially observed records.
Decision surface artifacts The single imputation methods introduce artifacts into the learnt decision surface. As the discriminative model tries to reflect the proper conditional probabilities p(y|xobs) at the location of the imputed values, the ambiguity about the class label needs to be reflected in the predictions. Notice how the MIWAE produces nearly optimal imputations in terms of minimizing a mean-square-error and how this translates into decision surface artifacts. These artifacts are related to the notion of imputation manifolds used by Le Morvan et al. (2021) to prove their main theorem on the Bayes consitency of impute-then-regress procedures. The remaining datasets are shown in figures 10–13 of appendix E, where also the permutation-invariance setup (appendix A.3) and histogram-based gradient boosting are included. Like the supMIWAE, these two approaches do not rely on explicit imputations.
As the supMIWAE can draw multiple importance samples, instead of one single imputation, the conditional probabilities p(y|xobs) can be properly reflected, when averaging over the importance samples. This avoids the single imputation artefacts.
Predictive performance Do the artifacts introduced by single imputation methods, compared to marginalizing over the missing values, translate into differences in predictive performance? The predictive performance corresponding to the decision surfaces in figure 3 are shown in table 1 in appendix E. Both the test-set accuracies and log likelihoods are more or less similar across the different strategies for handling missing values. This aligns with the results by Josse et al. (2019); Bertsimas et al. (2021); Le Morvan et al. (2021): the impute-then-regress procedures with powerful learners are consistent and Bayes optimal, i.e. when the learner has high capacity any single imputation can
essentially be undone. While single-imputation methods can lead to Bayes-optimal learners, Le Morvan et al. (2021) address that the difficulty of the learning problem depends on the choice of imputation function. This is seen in figure 4; while the supMIWAE and 0-imputation eventually ends up with similar performance the learning problem appears easier for the supMIWAE. For 0- imputation the learner has to introduce high-frequency artifacts to the decision surface in order to get the conditional probabilities p(y|xobs) calibrated, but in the case of a high capacity learner this can eventually be done.
Limited capacity learner We now turn to a situation where the capacity of the learner is limited in order to assess how single imputation methods fare against marginalizing over missing values. We start by noting that the decision surface for the Burger dataset, in the complete-data case, can be learnt using two sigmoids. When using MSE-optimal single imputations in the missing-values case this changes, as now 4 sigmoids are needed to make sure the imputations are well calibrated. This illustrates that the capacity of the learner determines whether single imputations can be undone. In figure 5 we explore how the capacity of the learner, in terms of the number of hidden units, affects the predictive performance, depending on the imputation strategy. More capacity experiments are found in figure 14.

heading: 4.2 IMAGE CLASSIFICATION
text: When the learning problems become harder, strong inductive biases are sometimes used to help learning. This is typically the case with image classification where convolutional neural networks (CNNs) are often used to introduce translational equivariance. The previous section showed that the strategy for handling missing data affects predictive performance when the learner is not in the high capacity regime. This is further investigated in terms of classification accuracy in image experiments with different MCAR missing processes, such as observed squares, missing squares and random dropout, see figures 15–17. The discriminator is a convolutional neural network, see appendix F.1, thus introducing an inductive bias in the learner that has proved useful in image classification in the complete-data case. While CNNs are fit for images, it should be noted that they are not universal approximators, and may consequently have a hard time undoing single imputations.
We apply the supMIWAE and single imputations methods to the MNIST dataset (LeCun et al., 1998) and the fashion MNIST dataset (Xiao et al., 2017) with three different MCAR missing mechanisms: 1) observed squares, where a 12×12 randomly placed square is observed, 2) missing squares, where a 12×12 randomly placed square is missing, and 3) dropout, where each pixel is missing with some constant probability m, over a range of missing rates.
In practice, we see that this setting where there is a strong inductive bias is similar to our previous low-capacity experiments: for both MNIST and Fashion MNIST, the supMIWAE is significantly
more accurate than single-imputation competitors, as seen in figure 6. Additional experiments, including experiments on natural images (SVHN) and visualizations of multiple imputations, are available in appendix F.

heading: 4.3 REGRESSION
text: We now turn to regression in smaller and lower dimensional datasets from the UCI database (Dua & Graff, 2017). An experimental setup modified from Hernández-Lobato & Adams (2015) and Skafte et al. (2019) is used here. Predictive performance in terms of root-mean-square-error is presented over a range of missing rates. Results are seen in figure 7, and more results can be found in figure 18 of appendix G. There is no clear winner here, although it appears that gradient boosting generally performs quite well, and that accurate single imputations (such as MIWAE’s or MissForest’s) generally outperform inaccurate ones (such as zero imputation). In general, our joint model performs as well as feeding accurate single imputations.

heading: 5 CONCLUSION
text: In the context of supervised deep learning with missing values, single imputation methods are an often used tool. While this approach is asymptotically consistent it incurs a bias in the learnt parameters. We have introduced a scalable approach to marginalizing over missing values using deep generative models, as a probabilistic alternative to single imputation. Besides limiting single imputation artifacts our results indicate that there are two cases where using a generative model can be quite valuable for supervised learning with missing values: when the classifier is not very flexible, or has a strong inductive bias. Moreover, leveraging such a generative model allows us to be able to use untouched any discriminative architecture (for instance, any pretrained one).

heading: A TRAINING DEEP SUPERVISED MODELS WITH MISSING DATA
text: In the following we group different strategies for handling missing data in supervised deep learning.
A.1 CONSTANT IMPUTATION
We denote imputation methods where there is a single constant for each feature which replaces missing data as constant imputation. We define an imputation function ι in general such that ι(x̃) ∈ X and ι(x̃)obs = xobs. A simple version of constant imputation is 0-imputation, which has the intuitive appeal that the activation from the input node is zeroed out (absent) and no gradient updates are made for the weights for that feature,
ι0(x̃) = x̃ s + 0 (1− s). (18)
In the unsupervised setting, constant imputation biases marginal and joint distributions, but Josse et al. (2019); Le Morvan et al. (2021); Bertsimas et al. (2021) have shown that mean imputation can be consistent in the supervised setting. Furthermore, Le Morvan et al. (2020b); Tresp et al. (1994) noted that the constants can be optimized with respect to the prediction loss.
Learnable parameters λ ∈ X can be inserted in place of the missing data, so that
ιλ(x̃) = x̃ s + λ (1− s). (19)
We denote this learnable imputation.
A.2 REGRESSIVE/MODEL BASED SINGLE IMPUTATION
In contrast to constant imputation, methods such as MICE (Buuren & Groothuis-Oudshoorn, 2010), MissForest (Stekhoven & Bühlmann, 2012), PPCA (Tipping & Bishop, 1999; Roweis, 1998) or flexible generative models (Mattei & Frellsen, 2019; Yoon et al., 2018) provides imputation of missing features conditional on observed features p(xmiss|xobs),
ι(x̃) = x̃ s + Ep(xmiss|xobs) [x] (1− s) (20)
In MICE, this conditional distribution comes from an implicit joint model over the covariates by regressing each covariate with missing values on all the others, in a round-robin fashion. MissForest is somewhat similar to MICE, using random forests as the regressors. Probabilistic PCA defines an explicit joint model over all covariates and single imputations of missing data can be obtained by conditioning on observed covariates. Here, MICE is a proper Bayesian multiple imputation method, while PPCA and MissForest draws multiple imputations without regard to parameter uncertainties. The recent advances in imputation of missing data using deep generative models can also be used here to the extent that they impose a joint distribution over covariates.
A.3 PERMUTATION INVARIANCE
Ma et al. (2018; 2019) used the permutation invariant setup (Zaheer et al., 2017; Qi et al., 2017) to handle missing data in the context of variational autoencoders. For each feature j ∈ obs(s), there is an embedding ej ∈ RU , which in general can hold any auxiliary information, but in our case will be a learnable embedding. Each element of the embedding is multiplied with the feature value tj = ej · xj and fed to a neural network h(·), mapping inputs from RU to RM where U is the embedding dimension and M is the dimension of the code-space being mapped into. This permutation invariant setup finally aggregates all the neural network outputs by a summation to get one fixed-length vector, that can be fed into a discriminative network,
ιPI(x obs) = ∑ j∈obs(s) h(tj). (21)
Note that in this case ιPI(xobs) ∈ RM , while in the other imputation functions ι(x̃) ∈ X .
A.4 TREE METHODS
Decision trees can in some cases handle missing data directly, without the need for imputation. There are different approaches to handle missing values in tree based methods, see (Josse et al., 2019) for a review. In our experiments we use histogram-based gradient boosting (Friedman, 2001), as implemented in scikit-learn (Pedregosa et al., 2011). This approach uses the Block Propagation method for the missing data. While a discriminative model using MLPs has a fixed capacity, gradient boosting methods can adapt its capacity in terms of tree depth based on a validation set error.

heading: B GRADIENTS
text: The log likelihood, when integrating over missing input features, is specified as
log pφ,θ(y,x obs) = log pφ(y|xobs) + log pθ(xobs) (22)
where log pφ(y|xobs) = log ∫ pφ(y|xobs,xmiss)pθ(xmiss|xobs) dxmiss. (23)
Both Tresp et al. (1994) and Ghahramani & Jordan (1995) analyzed the gradient of the log likelihood with respect to the discriminative network parameters,
∂ log pφ,θ(y,x obs) ∂φ = ∂ log pφ(y|xobs) ∂φ + ∂ log pθ(x obs) ∂φ (24)
= − 1 pφ(y|xobs) ∂pφ(y|xobs) ∂φ
(25)
= − 1 pφ(y|xobs) ∂ ∂φ
∫ pφ(y|xobs,xmiss)pθ(xmiss|xobs) dxmiss (26)
= − 1 pφ(y|xobs)
∫ ∂pφ(y|xobs,xmiss)
∂φ pθ(x
miss|xobs) dxmiss (27)
Here we see that the gradient of the discriminative network parameters, due to the classification/regression loss, is a weighted sum over all possible missing values, where the missing values are weighted according to their density in the data model, conditional on the observed data.
The gradient with respect to the θ parameters depends to some extend on the φ parameters.
∂ log pφ,θ(y,x obs) ∂θ = ∂ log pφ(y|xobs) ∂θ + ∂ log pθ(x obs) ∂θ (28)
= − 1 pφ(y|xobs) ∂pφ(y|xobs) ∂θ + ∂ log pθ(x obs) ∂θ (29)
= − 1 pφ(y|xobs) ∂ ∂θ
∫ pφ(y|xobs,xmiss)pθ(xmiss|xobs) dxmiss + ∂ log pθ(x obs)
∂θ (30)
= − 1 pφ(y|xobs) ∂ ∂θ Epθ(xmiss|xobs)
[ pφ(y|xobs,xmiss) ] + ∂ log pθ(x obs)
∂θ (31)
(32)
There are contributions to the gradient from both the generative part of the model and the discriminative part. This means that, if needed, the discriminative model can affect the generative model, as is the case in (Ipsen et al., 2021). The gradient of the expectation can be approached using the reparameterization trick (Kingma & Welling, 2014; Rezende et al., 2014).

heading: C RELATION TO MULTIPLE IMPUTATION
text: The gold standard for imputation is multiple imputation (Rubin, 1996). This requires learning a model of the observed covariates xobs with parameters θ. K sets of imputations are then drawn from
the posterior predictive distribution
xmissk ∼ p(xmiss|xobs), k = 1, . . . ,K (33) where p(xmiss|xobs) = ∫ p(xmiss|θ)p(θ|xobs) dθ if the model is Bayesian (see e.g. Gelman et al., 2013, chapter 18.2). Several frequentist multiple imputation methods also exist, where the uncertainty with respect to the imputation model parameters is disregarded. In that case, the Bayesian posterior predictive is replaced by a conditional distribution p(xmiss|xobs) = pθ̂(x
miss|xobs), where θ̂ is any point estimate of θ.
In the case of regression/classification, if we roughly follow the notations of our paper, the analysis proceeds as follows
• For each set of imputed values (xmissik )i≤n,k≤K , use a learning algorithm to estimate the discriminative parameters φ1, . . . ,φK :
φ̂1 = arg max φ n∑ i=1 log pφ(yi|xobsi ,xmissi1 ) (34)
...
φ̂K = arg max φ n∑ i=1 log pφ(yi|xobsi ,xmissiK ) (35)
• Obtain an average estimate of φ using
φ̄K = 1
K K∑ k=1 φ̂k (36)
and a variance estimate can similarly be obtained as described in (Gelman et al., 2013; Rubin, 1996).
In the limit of infinite imputations, we can write this as the conditional expectation
φ̄ = Exmiss1 ,...,xmissn arg max φ n∑ i=1 log pφ(yi|xobsi ,xmissi ) ∣∣∣∣∣∣xobs1 , . . . ,xobsn  . (37)
Instead of fitting K conditional models to K imputed datasets, our approach (when θ is fixed) is, in comparison, to directly maximize the expectation of the conditional distribution with respect to the missing data
φ̄ = arg max φ n∑ i=1 logExmissi [ pφ(yi|xobsi ,xmissi )|xobsi ] . (38)
Also note that we have no need for several fixed sets of imputations, at each gradient step a new set of samples are drawn from the imputation model, xmiss ∼ pθ(xmiss|xobs). The number of imputations can be set dynamically, using ∼ 5–50 during training and ∼ 10k for testing. Our approach, using DLVMs as imputation models, is not Bayesian, but any available proper Bayesian imputation model could be used in place of the DLVMs.
In practice, an important difference between traditional multiple imputation and our approach, is that multiple imputation uses a uniform weight 1/K for each imputation, while we use importance weights, that will treat each imputation based on its likelihood.

heading: D THEORETICAL PROPERTIES OF THE SUPMIWAE BOUND
text: We will briefly see how the supMIWAE bound has the same theoretical properties as IWAE bounds. We follow Ipsen et al. (2021, appendix D), who were following Burda et al. (2016) and Domke & Sheldon (2018). Recall the definition of the bound: we have
RK(x obs,y) =
1
K K∑ k=1 pφ(y|xobs,xmissk )pθ(xobs|zk)pθ(zk) qγ(zk|xobs) , (39)
and
LK(θ,φ,γ) = n∑ i=1 E(xmissk ,zk) [ logRK(x obs i ,yi) ] . (40)
The ith term of the sum can be seen as a regular IWAE bound with observation (xobs,y), latent variable (z,xmiss) with prior pθ(xmiss|z)pθ(z), and variational distribution pθ(xmiss|z)qγ(z|xobs). Therefore, Theorem 1 from Burda et al. (2016) can be applied, and leads to the monotonicity of the bound: L1(θ,φ,γ) ≤ . . . ≤ LK(θ, φ, γ). To show that the bound additionally converges to the log-likelihood when K −→ ∞, we use Theorem 3 of Domke & Sheldon (2018) for all n terms of the sum that constitutes the bound. Theorem. Assuming that, for all i ∈ {1, ..., n},
• there exists αi > 0 such that E [ |R1(xobsi ,yi)− pθ,φ(xobsi ,yi)|2+αi ] <∞,
• lim supK−→∞ E [ 1/RK ] <∞,
the supMIWAE bound converges to the true joint likelihood at rate 1/K: n∑ i=1 log pφ,θ(yi,x obs i )− LK(θ,φ,γ) ∼ K→∞ 1 K n∑ i=1 Var[R1(xobsi ,yi)] 2pθ,φ(x obs i ,yi) 2 . (41)
E 2D CLASSIFICATION
The datasets used in the 2D classification experiments are shown in figure 8. The Half-moons and Circles datasets are from Scikit-learn (Pedregosa et al., 2011). The Pin-wheel dataset is modified from Johnson et al. (2016)1. The Burger dataset. we have designed to illustrate the case where optimal single imputation requires severe changes in the decision surface.
The discriminative model is an MLP with 3 hidden layers with 50 hidden units each. The final layer parameterizes a categorical distribution over the dataset classes. The generative model used for the MIWAE single imputations and the supMIWAE consists of an encoder and decoder with 3 hidden layers containing 50 hidden units each. The encoder parameterizes an isotropic Gaussian distribution over a latent space with dimension p, same as the feature space. The deocder parameterizes an isotropic Gaussian distribution over the feature space. In figure 9 samples from the generative model are shown, that is, samples from the MIWAE trained on each of the datasets. In figures 10–13 decision surfaces and imputations are shown. The corresponding quantitative results are in tables 1– 4.
1https://github.com/mattjj/svae
E.1 CAPACITY
The predictive performance of the discriminator at different capacities, in terms of number of hidden units, is shown in figure 14.
F IMAGE CLASSIFICATION
F.1 TRAINING DETAILS FOR MNIST/FMNIST
For all MNIST/fMNIST experiments, the discriminative model is a convolutional neural network with four layers and a categorical distribution over classes, see table 5. The generative model uses a convolutional architecture similar to the DCGAN (Radford et al., 2015) and a latent space of dimension 128. A Continuous Bernoulli (Loaiza-Ganem & Cunningham, 2019) observation model is used for MNIST and the discretized logistic distribution (Salimans et al., 2017) is used as the observation model for Fashion MNIST. During training, K = 5 importance samples and a batch size of 128 are used.
F.2 PRELIMINARY RESULTS ON NATURAL IMAGES (SVHN)
In order to assess classification of natural images we now turn to the Street View House Numbers dataset (SVHN, Netzer et al., 2011), with the “observed squares” missing mechanism, where 20×20 randomly placed squares are observed. We apply both a low capacity (4 layer) discriminator and a higher capacity (8 layers) using the convolutional gated blocks from Dauphin et al. (2017), see architectures in tables 7 and 8. In table 6, test set accuracies are shown and for the 8 layer discriminator, the accuracy on the fully observed test set is also shown. In the low capacity regime, the supMIWAE outperforms the two single imputation methods and performs on par with MIWAE. With the higher capacity discriminator supMIWAE has performance similar to 0-imputation and learnable imputation, while the MIWAE imputation is slightly worse. This displays the same phenomenon as in the 2D classification experiments, where powerful learners are able to undo single imputations. While the discriminative models learnt based on constant imputation can achieve the same performance as the supMIWAE in the missing case, the kinds of discriminative models that are learnt are different. When applied to a fully observed test-set the discriminator learnt by the supMIWAE outperforms the other models.

heading: G REGRESSION
text: All discriminative models are neural networks with one hidden layer and 50 hidden units, where the final layer of the neural networks parameterizes a Gaussian distribution. The datasets are split randomly 20 times with 90% of the data in a training set and 10% in a test set. A validation set with 10% of the training data is used for early stopping and a batch size of 256 is used. For the permutation invariant setup, an embedding of dimension U = 20 and fixed length vectors ιPI(xobs) of dimensionM = 10 are used. Missing values are introduced completely at random in the training, validation and test sets by removing each feature with probability m, where m is the missing rate. In order to avoid reducing the sample size, due to completely missing covariates, whenever all the covariates of an observation go missing, we set one of them to observed, selected at random.
For the generative models, 2 layered neural networks with 128 hidden units are used for the encoder and the decoder and a latent space of dimension 10. A Gaussian observation model is used and K = 25 importance samples are used during training.
We compare different imputation techniques, modelling approaches and a non-deep competitor. Zero-imputation, MissForest (Stekhoven & Bühlmann, 2012) and MICE (Buuren & GroothuisOudshoorn, 2010) are different imputation techniques applied, before training the discriminative model. Learnable imputation is using learnable constants as imputation, while permutation invariance is utilizing a learnable embedding. Finally, histogram-based gradient boosting (Friedman, 2001), as implemented in sklearn (Pedregosa et al., 2011), is a non-deep baseline, using the validation set for early stopping. Results are seen in figure 18.

Review summaries:
The paper handles the issue of missing values in supervised deep learning settings. The fig.1 describes their method aptly. Their method (supMIWAE) is a combination of a VAE with a neural network classifier. Given the task to predict p(Y|x_{obs}, x_{miss}), the authors view it as a joint model of covariates and outcomes. Their model takes in x_{obs} and first fits a distribution to get x_{miss} and then models them as a joint distribution to predict outcome Y|(x_{obs}, x_{miss}). They use importance sampling technique to get multiple samples from the generative model. 

This paper approaches the problem of supervised learning with missing data. The authors propose a probabilistic approach by jointly modeling the observed data, missing data and outcomes. The main contribution of this model is that they rely on deep generative models which seems to be an improvement from previous methods based on simpler generative models like mixture of Gaussians, for example. With the proposed model, the authors derive an optimization problem that consists in optimizing the discriminative model (classifier) and the generative model simultaneously. For the training and testing phases, the proposed algorithm considers multiple imputations of missing data, which is demonstrated to be superior to single imputation methods. Experimental results on small and simple datasets (2D datasets, MNIST digits-Fashion and regression) are nicely presented and analyzed.