Year: 2021

Abstract:
When a missing process depends on the missing values themselves, it needs to be explicitly modelled and taken into account while doing likelihood-based inference. We present an approach for building and fitting deep latent variable models (DLVMs) in cases where the missing process is dependent on the missing data. Specifically, a deep neural network enables us to flexibly model the conditional distribution of the missingness pattern given the data. This allows for incorporating prior information about the type of missingness (e.g. self-censoring) into the model. Our inference technique, based on importance-weighted variational inference, involves maximising a lower bound of the joint likelihood. Stochastic gradients of the bound are obtained by using the reparameterisation trick both in latent space and data space. We show on various kinds of data sets and missingness patterns that explicitly modelling the missing process can be invaluable.

text: 1 INTRODUCTION
z
x
s
θ
φ
γ
N
(a)
PPCA not-MIWAE PPCA
(b)
Figure 1: (a) Graphical model of the not-MIWAE. (b) Gaussian data with MNAR values. Dots are fully observed, partially observed data are displayed as black crosses. A contour of the true distribution is shown together with directions found by PPCA and not-MIWAE with a PPCA decoder.
Missing data often constitute systemic issues in real-world data analysis, and can be an integral part of some fields, e.g. recommender systems. This requires the analyst to take action by either using methods and models that are applicable to incomplete data or by performing imputations of the missing data before applying models requiring complete data. The expected model performance (often measured in terms of imputation error or innocuity of missingness on the inference results) depends on the assumptions made about the missing mechanism and how well those assumptions match the true missing mechanism. In a seminal paper, Rubin (1976) introduced a formal probabilistic framework to assess missing mechanism assumptions and their consequences. The most commonly used assumption, either implicitly or explicitly, is that a part of the data is missing at random
(MAR). Essentially, the MAR assumption means that the missing pattern does not depend on the missing values. This makes it possible to ignore the missing data mechanism in likelihood-based inference by marginalizing over the missing data. The often implicit assumption made in nonprobabilistic models and ad-hoc methods is that the data are missing completely at random (MCAR). MCAR is a stronger assumption than MAR, and informally it means that both observed and missing data do not depend on the missing pattern. More details on these assumptions can be found in the monograph of Little & Rubin (2002); of particular interest are also the recent revisits of Seaman et al. (2013) and Doretti et al. (2018). In this paper, our goal is to posit statistical models that leverage deep learning in order to break away from these assumptions. Specifically, we propose a general ∗Department of Applied Mathematics and Computer Science, Technical University of Denmark, Denmark †Université Côte d’Azur, Inria (Maasai team), Laboratoire J.A. Dieudonné, UMR CNRS 7351, France ‡Equal contribution
recipe for dealing with cases where there is prior information about the distribution of the missing pattern given the data (e.g. self-censoring).
The MAR and MCAR assumptions are violated when the missing data mechanism is dependent on the missing data themselves. This setting is called missing not at random (MNAR). Here the missing mechanism cannot be ignored, doing so will lead to biased parameter estimates. This setting generally requires a joint model for data and missing mechanism.
Deep latent variable models (DLVMs, Kingma & Welling, 2013; Rezende et al., 2014) have recently been used for inference and imputation in missing data problems (Nazabal et al., 2020; Ma et al., 2018; 2019; Ivanov et al., 2019; Mattei & Frellsen, 2019). This led to impressive empirical results in the MAR and MCAR case, in particular for high-dimensional data.

heading: 1.1 CONTRIBUTIONS
text: We introduce the not-missing-at-random importance-weighted autoencoder (not-MIWAE) which allows for the application of DLVMs to missing data problems where the missing mechanism is MNAR. This is inspired by the missing data importance-weighted autoencoder (MIWAE, Mattei & Frellsen, 2019), a framework to train DLVMs in MAR scenarios, based itself on the importanceweighted autoencoder (IWAE) of Burda et al. (2016). The general graphical model for the notMIWAE is shown in figure 1a. The first part of the model is simply a latent variable model: there is a stochastic mapping parameterized by θ from a latent variable z ∼ p(z) to the data x ∼ pθ(x|z), and the data may be partially observed. The second part of the model, which we call the missing model, is a stochastic mapping from the data to the missing mask s ∼ pφ(s|x). Explicit specification of the missing model pφ(s|x) makes it possible to address MNAR issues. The model can be trained efficiently by maximising a lower bound of the joint likelihood (of the observed features and missing pattern) obtained via importance weighted variational inference (Burda et al., 2016). A key difference with the MIWAE is that we use the reparameterization trick in the data space, as well as in the code space, in order to get stochastic gradients of the lower bound.
Missing processes affect data analysis in a wide range of domains and often the MAR assumption does not hold. We apply our method to censoring in datasets from the UCI database, clipping in images and the issue of selection bias in recommender systems.

heading: 2 BACKGROUND
text: Assume that the complete data are stored within a data matrix X = (x1, . . . ,xn)ᵀ ∈ Xn that contain n i.i.d. copies of the random variable x ∈ X , where X = X1 × · · · × Xp is a p-dimensional feature space. For simplicity, xij refers to the j’th feature of xi, and xi refers to the i’th sample in the data matrix. Throughout the text, we will make statements about the random variable x, and only consider samples xi when necessary. In a missing data context, each sample can be split into an observed part and a missing part, xi = (xoi ,x m i ). The pattern of missingness is individual to each copy of x and described by a corresponding mask random variable s ∈ {0, 1}p. This leads to a mask matrix S = (s1, . . . , sn)ᵀ ∈ {0, 1}n×p verifying sij = 1 if xij is observed and sij = 0 if xij is missing.
We wish to construct a parametric model pθ,φ(x, s) for the joint distribution of a single data point x and its mask s, which can be factored as
pθ,φ(x, s) = pθ(x)pφ(s|x). (1)
Here pφ(s|x) = pφ(s|xo,xm) is the conditional distribution of the mask, which may depend on both the observed and missing data, through its own parameters φ. The three assumptions from the framework of Little & Rubin (2002) (see also Ghahramani & Jordan, 1995) pertain to the specific form of this conditional distribution:
• MCAR: pφ(s|x) = pφ(s), • MAR: pφ(s|x) = pφ(s|xo), • MNAR: pφ(s|x) may depend on both xo and xm.
To maximize the likelihood of the parameters (θ, φ), based only on observed quantities, the missing data is integrated out from the joint distribution
pθ,φ(x o, s) = ∫ pθ(x o,xm)pφ(s|xo,xm) dxm. (2)
In both the MCAR and MAR cases, inference for θ using the full likelihood becomes proportional to pθ,φ(xo, s) ∝ pθ(xo), and the missing mechanism can be ignored while focusing only on pθ(xo). In the MNAR case, the missing mechanism can depend on both observed and missing data, offering no factorization of the likelihood in equation (2). The parameters of the data generating process and the parameters of the missing data mechanism are tied together by the missing data.

heading: 2.1 PPCA EXAMPLE
text: A linear DLVM with isotropic noise variance can be used to recover a model similar to probabilistic principal component analysis (PPCA, Roweis, 1998; Tipping & Bishop, 1999). In figure 1b, a dataset affected by an MNAR missing process is shown together with two fitted PPCA models, regular PPCA and the not-MIWAE formulated as a PPCA-like model. Data is generated from a multivariate normal distribution and an MNAR missing process is imposed by setting the horizontal coordinate to missing when it is larger than its mean, i.e. it becomes missing because of the value it would have had, had it been observed. Regular PPCA for missing data assumes that the missing mechanism is MAR so that the missing process is ignorable. This introduces a bias, both in the estimated mean and in the estimated principal signal direction of the data. The not-MIWAE PPCA assumes the missing mechanism is MNAR so the data generating process and missing data mechanism are modelled jointly as described in equation (2).

heading: 2.2 PREVIOUS WORK
text: In (Rubin, 1976) the appropriateness of ignoring the missing process when doing likelihood based or Bayesian inference was introduced and formalized. The introduction of the EM algorithm (Dempster et al., 1977) made it feasible to obtain maximum likelihood estimates in many missing data settings, see e.g. Ghahramani & Jordan (1994; 1995); Little & Rubin (2002). Sampling methods such as Markov chain Monte Carlo have made it possible to sample a target posterior in Bayesian models, including the missing data, so that parameter marginal distributions and missing data marginal distributions are available directly (Gelman et al., 2013). This is also the starting point of the multiple imputations framework of Rubin (1977; 1996). Here the samples of the missing data are used to provide several realisations of complete datasets where complete-data methods can be applied to get combined mean and variability estimates.
The framework of Little & Rubin (2002) is instructive in how to handle MNAR problems and a recent review of MNAR methods can be found in (Tang & Ju, 2018). Low rank models were used for estimation and imputation in MNAR settings by Sportisse et al. (2020a). Two approaches were taken to fitting models, 1) maximising the joint distribution of data and missing mask using an EM algorithm, and 2) implicitly modelling the joint distribution by concatenating the data matrix and the missing mask and working with this new matrix. This implies a latent representation both giving rise to the data and the mask. An overview of estimation methods for PCA and PPCA with missing data was given by Ilin & Raiko (2010), while PPCA in the presence of an MNAR missing mechanism has been addressed by Sportisse et al. (2020b). There has been some focus on MNAR issues in the form of selection bias within the recommender system community (Marlin et al., 2007; Marlin & Zemel, 2009; Steck, 2013; Hernández-Lobato et al., 2014; Schnabel et al., 2016; Wang et al., 2019) where methods applied range from joint modelling of data and missing model using multinomial mixtures and matrix factorization to debiasing existing methods using propensity based techniques from causality.
Deep latent variable models are intuitively appealing in a missing context: the generative part of the model can be used to sample the missing part of an observation. This was already utilized by Rezende et al. (2014) to do imputation and denoising by sampling from a Markov chain whose stationary distribution is approximately the conditional distribution of the missing data given the observed. This procedure has been enhanced by Mattei & Frellsen (2018a) using Metropolis-withinGibbs. In both cases the experiments were assuming MAR and a fitted model, based on complete data, was already available.
Approaches to fitting DLVMs in the presence of missing have recently been suggested, such as the HI-VAE by Nazabal et al. (2020) using an extension of the variational autoencoder (VAE) lower bound, the p-VAE by Ma et al. (2018; 2019) using the VAE lower bound and a permutation invariant encoder, the MIWAE by Mattei & Frellsen (2019), extending the IWAE lower bound (Burda et al., 2016), and GAIN (Yoon et al., 2018) using GANs for missing data imputation. All approaches are assuming that the missing process is MAR or MCAR. In (Gong et al., 2020), the data and missing mask are modelled together, as both being generated by a mapping from the same latent space, thereby tying the data model and missing process together. This gives more flexibility in terms of missing process assumptions, akin to the matrix factorization approach by Sportisse et al. (2020a).
In concurrent work, Collier et al. (2020) have developed a deep generative model of the observed data conditioned on the mask random variable, and Lim et al. (2021) apply a model similar to the not-MIWAE to electronic health records data. In forthcoming work, Ghalebikesabi et al. (2021) propose a deep generative model for non-ignorable missingness building on ideas from VAEs and pattern-set mixture models.

heading: 3 INFERENCE IN DLVMS AFFECTED BY MNAR
text: In an MNAR setting, the parameters for the data generating process and the missing data mechanism need to be optimized jointly using all observed quantities. The relevant quantity to maximize is therefore the log-(joint) likelihood
`(θ, φ) = n∑ i=1 log pθ,φ(x o i , si), (3)
where we can rewrite the general contribution of data points log pθ,φ(xo, s) as
log ∫ pφ(s|xo,xm)pθ(xo|z)pθ(xm|z)p(z) dz dxm, (4)
using the assumption that the observation model is fully factorized pθ(x|z) = ∏ j pθ(xj |z), which implies pθ(x|z) = p(xo|z)pθ(xm|z). The integrals over missing and latent variables make direct maximum likelihood intractable. However, the approach of Burda et al. (2016), using an inference network and importance sampling to derive a more tractable lower bound of `(θ, φ), can be used here as well. The key idea is to posit a conditional distribution qγ(z|xo) called the variational distribution that will play the role of a learnable proposal in an importance sampling scheme.
As in VAEs (Kingma & Welling, 2013; Rezende et al., 2014) and IWAEs (Burda et al., 2016), the distribution qγ(z|xo) comes from a simple family (e.g. the Gaussian or Student’s t family) and its parameters are given by the output of a neural network (called inference network or encoder) that takes xo as input. The issue is that a neural net cannot readily deal with variable length inputs (which is the case of xo). This was tackled by several works: Nazabal et al. (2020) and Mattei & Frellsen (2019) advocated simply zero-imputing xo to get inputs with constant length, and Ma et al. (2018; 2019) used a permutation-invariant network able to deal with inputs with variable length.
Introducing the variational distribution, the contribution of a single observation is equal to
log pθ,φ(x o, s) = log
∫ pφ(s|xo,xm)pθ(xo|z)p(z)
qγ(z|xo) qγ(z|xo)pθ(xm|z) dxm dz (5)
= logEz∼qγ(z|xo),xm∼pθ(xm|z)
[ pφ(s|xo,xm)pθ(xo|z)p(z)
qγ(z|xo)
] . (6)
The main idea of importance weighed variational inference and of the IWAE is to replace the expectation inside the logarithm by a Monte Carlo estimate of it (Burda et al., 2016). This leads to the objective function
LK(θ, φ, γ) = n∑ i=1 E log 1 K K∑ k=1 wki  , (7) where, for all k ≤ K, i ≤ n,
wki = pφ(si|xoi ,xmki)pθ(xoi |zki)p(zki)
qγ(zki|xoi ) , (8)
and (z1i,xm1i), . . . , (zKi,x m Ki) are K i.i.d. samples from qγ(z|xoi )pθ(xm|z), over which the expectation in equation (7) is taken. The unbiasedness of the Monte Carlo estimates ensures (via Jensen’s inequality) that the objective is indeed a lower-bound of the likelihood. Actually, under the moment conditions of (Domke & Sheldon, 2018, Theorem 3), which we detail in Appendix D, it is possible to show that the sequence (LK(θ, φ, γ))K≥1 converges monotonically (Burda et al., 2016, Theorem 1) to the likelihood:
L1(θ, φ, γ) ≤ . . . ≤ LK(θ, φ, γ) −−−−→ K→∞ `(θ, φ). (9)
Properties of the not-MIWAE objective The bound LK(θ, φ, γ) has essentially the same properties as the (M)IWAE bounds, see Mattei & Frellsen, 2019, Section 2.4 for more details. The key difference is that we are integrating over both the latent space and part of the data space. This means that, to obtain unbiased estimates of gradients of the bound, we will need to backpropagate through samples from qγ(z|xoi )pθ(xm|z). A simple way to do this is to use the reparameterization trick both for qγ(z|xoi ) and pθ(xm|z). This is the approach that we chose in our experiments. The main limitation is that pθ(x|z) has to belong to a reparameterizable family, like Gaussians or Student’s t distributions (see Figurnov et al., 2018 for a list of available distributions). If the distribution is not readily reparametrisable (e.g. if the data are discrete), several other options are available, see e.g. the review of Mohamed et al. (2020), and, in the discrete case, the continuous relaxations of Jang et al. (2017) and Maddison et al. (2017).
Imputation When the model has been trained, it can be used to impute missing values. If our performance metric is a loss function L(xm, x̂m), optimal imputations x̂m minimise Exm [L(xm, x̂m)|xo, s]. When L is the squared error, the optimal imputation is the conditional mean that can be estimated via self-normalised importance sampling (Mattei & Frellsen, 2019), see appendix B for more details.

heading: 3.1 USING PRIOR INFORMATION VIA THE MISSING DATA MODEL
text: The missing data mechanism can both be known/decided upon in advance (so that the full relationship pφ(s|x) is fixed and no parameters need to be learned) or the type of missing mechanism can be known (but the parameters need to be learnt) or it can be unknown both in terms of parameters and model. The more we know about the nature of the missing mechanism, the more information we can put into designing the missing model. This in turn helps inform the data model how its parameters should be modified so as to accommodate the missing model. This is in line with the findings of Molenberghs et al. (2008), who showed that, for MNAR modelling to work, one has to leverage prior knowledge about the missing process. A crucial issue is under what model assumptions the full data distribution can be recovered from incomplete sample. Indeed, some general missing models may lead to inconsistent statistical estimation (see e.g. Mohan & Pearl, 2021; Nabi et al., 2020).
The missing model is essentially solving a classification problem; based on the observed data and the output from the data model filling in the missing data, it needs to improve its “accuracy” in predicting the mask. A Bernoulli distribution is used for the probability of the mask given both observed and missing data
pφ(s|xo,xm) = pφ(s|x) = Bern(s|πφ(x)) = ∏p j=1 πφ,j(x) sj (1− πφ,j(x))1−sj . (10)
Here πj is the estimated probability of being observed for that particular observation for feature j. The mapping πφ,j(x) from the data to the probability of being observed for the j’th feature can be as general or specific as needed. A simple example could be that of self-masking or self-censoring, where the probability of the j’th feature being observed is only dependent on the feature value, xj . Here the mapping can be a sigmoid on a linear mapping of the feature value, πφ,j(x) = σ(axj + b). The missing model can also be based on a group theoretic approach, see appendix C.

heading: 4 EXPERIMENTS
text: In this section we apply the not-MIWAE to problems with values MNAR: censoring in multivariate datasets, clipping in images and selection bias in recommender systems. Implementation details and a link to source code can be found in appendix A.

heading: 4.1 EVALUATION METRICS
text: Model performance can be assessed using different metrics. A first metric would be to look at how well the marginal distribution of the data has been inferred. This can be assessed, if we happen to have a fully observed test-set available. Indeed, we can look at the test log-likelihood of this fully observed test-set as a measure of how close pθ(x) and the true distribution of x are. In the case of a DLVM, performance can be estimated using importance sampling with the variational distribution as proposal (Rezende et al., 2014). Since the encoder is tuned to observations with missing data, it should be retrained (while keeping the decoder fixed) as suggested by Mattei & Frellsen (2018b).
Another metric of interest is the imputation error. In experimental settings where the missing mechanism is under our control, we have access to the actual values of the missing data and the imputation error can be found directly as an error measure between these and the reconstructions from the model. In real-world datasets affected by MNAR processes, we cannot use the usual approach of doing a train-test split of the observed data. As the test-set is biased by the same missing mechanism as the training-set it is not representative of the full population. Here we need a MAR data sample to evaluate model performance (Marlin et al., 2007).

heading: 4.2 SINGLE IMPUTATION IN UCI DATA SETS AFFECTED BY MNAR
text: We compare different imputation techniques on datasets from the UCI database (Dua & Graff, 2017), where in an MCAR setting the MIWAE has shown state of the art performance (Mattei & Frellsen, 2019). An MNAR missing process is introduced by self-masking in half of the features: when the feature value is higher than the feature mean it is set to missing. The MIWAE and not-MIWAE, as well as their linear PPCA-like versions, are fitted to the data with missing values. For the notMIWAE three different approaches to the missing model are used: 1) agnostic where the data model output is mapped to logits for the missing process via a single dense linear layer, 2) self-masking where logistic regression is used for each feature and 3) self-masking known where the sign of the weights in the logistic regression is known.
We compare to the low-rank approximation of the concatenation of data and mask by Sportisse et al. (2020a) that is implicitly modelling the data and mask jointly. Furthermore we compare to mean imputation, missForest (Stekhoven & Bühlmann, 2012) and MICE (Buuren & GroothuisOudshoorn, 2010) using Bayesian Ridge regression. Similar settings are used for the MIWAE and not-MIWAE, see appendix A. Results over 5 runs are seen in table 1. Results for varying missing rates are in appendix E.
The low-rank joint model is almost always better than PPCA, missForest, MICE and mean, i.e. all M(C)AR approaches, which can be attributed to the implicit modelling of data and mask together. At the same time the not-MIWAE PPCA is always better than the corresponding low-rank joint model, except for the agnostic missing model on the Yeast dataset. Supplying the missing model with more knowledge of the missing process (that it is self-masking and the direction of the missing mechanism) improves performance. The not-MIWAE performance is also improved with more knowledge in the missing model. The agnostic missing process can give good performance, but is
often led astray by an incorrectly learned missing model. This speaks to the trade-off between data model flexibility and missing model flexibility. The not-MIWAE PPCA has huge inductive bias in the data model and so we can employ a more flexible missing model and still get good results. For the not-MIWAE having both a flexible data model and a flexible missing model can be detrimental to performance. One way to asses the learnt missing processes is the mask classification accuracy on fully observed data. These are reported in table A1 and show that the accuracy increases as more information is put into the missing model.

heading: 4.3 CLIPPING IN SVHN IMAGES
text: We emulate the clipping phenomenon in images on the street view house numbers dataset (SVHN, Netzer et al., 2011). Here we introduce a self-masking missing mechanism that is identical for all pixels. The missing data is Bernoulli sampled with probability
Pr(sij = 1|xij) = 1
1 + e−logits , logits =W (xij − b), (11)
where W = −50 and b = 0.75. This mimmicks a clipping process where 0.75 is the clipping point (the data is converted to gray scale in the [0, 1] range). For this experiment we use the true missing process as the missing model in the not-MIWAE.
Table 2 shows model performance in terms of imputation RMSE and test-set log likelihood as estimated with 10k importance samples. The not-MIWAE outperforms the MIWAE both in terms of test-set log likelihood and imputation RMSE. This is further illustrated in the imputations shown in figure 3. Since the MIWAE is only fitting the observed data, the range of pixel values in the imputations is limited compared to the true range. The not-MIWAE is forced to push some of the data-distribution towards higher pixel values, in order to get a higher likelihood in the logistic regression in the missing model. In figures 2a–2c, histograms over the imputation values are shown together with the true pixel values of the missing data. Here we see that the not-MIWAE puts a considerable amount of probability mass above the clipping value.

heading: 4.4 SELECTION BIAS IN THE YAHOO! R3 DATASET
text: The Yahoo! R3 dataset (webscope.sandbox.yahoo.com) contains ratings on a scale from 1–5 of songs in the database of the Yahoo! LaunchCast internet radio service and was first presented in (Marlin et al., 2007). It consists of two datasets with the same 1,000 songs selected randomly from
the LaunchCast database. The first dataset is considered an MNAR training set and contains selfselected ratings from 15,400 users. In the second dataset, considered an MCAR test-set, 5,400 of these users were asked to rate exactly 10 randomly selected songs. This gives a unique opportunity to train a model on a real-world MNAR-affected dataset while being able to get an unbiased estimate of the imputation error, due to the availability of MCAR ratings. The plausibility that the set of selfselected ratings was subject to an MNAR missing process was explored and substantiated by Marlin et al. (2007). The marginal distributions of samples from the self-selected dataset and the randomly selected dataset can be seen in figures 4a and 4b.
We train the MIWAE and the not-MIWAE on the MNAR ratings and evaluate the imputation error on the MCAR ratings. Both a gaussian and a categorical observation model is explored. In order to get reparameterized samples in the data space for the categorical observation model, we use the Gumbel-Softmax trick (Jang et al., 2017) with a temperature of 0.5. The missing model is a logistic regression for each item/feature, with a shared weight across features and individual biases. A description of competitors can be found in appendix A.3 and follows the setup in (Wang et al., 2019). The results are grouped in table 3, from top to bottom, according to models not including the missing process (MAR approaches), models using propensity scoring techniques to debias training losses, and finally models learning a data model and a missing model jointly, without the use of propensity estimates.
The not-MIWAE shows state of the art performance, also compared to models based on propensity scores. The propensity based techniques need access to a small sample of MCAR data, i.e. a part of the test-set, to estimate the propensities using Naive Bayes, though they can be estimated using logistic regression if covariates are available (Schnabel et al., 2016) or using a nuclear-norm-constrained matrix factorization of the missing mask itself (Ma & Chen, 2019). We stress that the not-MIWAE does not need access to similar unbiased data in order to learn the missing model. However, the missing model in the not-MIWAE can take available information into account, e.g. we could fit a continuous mapping to the propensities and use this as the missing model, if propensities were available. Histograms over imputations for the missing data in the MCAR test-set can be seen for the MIWAE and notMIWAE in figures 4c and 4d. The marginal distribution of the not-MIWAE imputations are seen to match that of the MCAR test-set better than the marginal distribution of the MIWAE imputations.

heading: 5 CONCLUSION
text: The proposed not-MIWAE is versatile both in terms of defining missing mechanisms and in terms of application area. There is a trade-off between data model complexity and missing model complexity. In a parsimonious data model a very general missing process can be used while in flexible data
model the missing model needs to be more informative. Specifically, any knowledge about the missing process should be incorporated in the missing model to improve model performance. Doing so using recent advances in equivariant/invariant neural networks is an interesting avenue for future research (see appendix C). Recent developments on the subject of recoverability/identifiability of MNAR models (Sadinle & Reiter, 2018; Mohan & Pearl, 2021; Nabi et al., 2020; Sportisse et al., 2020b) could also be leveraged to design provably idenfiable not-MIWAE models.
Several extensions of the graphical models used here could be explored. For example, one could break off the conditional independence assumptions, in particular the one of the mask given the data. This could, for example, be done by using an additional latent variable pointing directly to the mask. Combined with a discriminative classifier, the not-MIWAE model could also be used in supervised learning with input values missing not at random following the techniques by Ipsen et al. (2020).

heading: ACKNOWLEDGMENTS
text: The Danish Innovation Foundation supported this work through Danish Center for Big Data Analytics driven Innovation (DABAI). JF acknowledge funding from the Independent Research Fund Denmark (grant number 9131-00082B) and the Novo Nordisk Foundation (grant numbers NNF20OC0062606 and NNF20OC0065611).

heading: C MISSING MODEL, GROUP THEORETIC APPROACH
text: A more complex form of prior information that can be used to choose the form of πφ(x) is grouptheoretic. For example, we may know a priori that pφ(s|x) is invariant to a certain group action g ·x on the data space:
∀g, pφ(s|x) = pφ(s|g · x). (20)
This would for example be the case, if the data sets were made of images whose class is invariant to translations (which is the case of most image data sets, like MNIST or SVHN), and with a missing model only dependent on the class. Similarly, one may know that the missing process is equivariant:
∀g, pφ(g · s|x) = pφ(s|g−1 · x). (21)
Again, such a setting can appear when there is strong geometric structure in the data (e.g. with images or proteins). Invariance or equivariance can be built in the architecture of πφ(x) by leveraging the quite large body of work on invariant/equivariant convolutional neural networks, see e.g. Bietti & Mairal (2017); Cohen et al. (2019); Zaheer et al. (2017); Wiqvist et al. (2019); Bloem-Reddy & Teh (2020), and references therein.

heading: D THEORETICAL PROPERTIES OF THE NOT-MIWAE BOUND
text: The properties of the not-MIWAE bound are directly inherited from the ones of the usual IWAE bound. Indeed, as we will see, the not-MIWAE bound is a particular instance of IWAE bound with an extended latent space composed of both the code and the missing values. More specifically, recall the definition of the not-MIWAE bound
LK(θ, φ, γ) = n∑ i=1 E log 1 K K∑ k=1 wki  , with wki = pθ(xoi |zki)pφ(si|xoi ,xmki)p(zki) qγ(zki|xoi ) . (22)
Each ith term of the sum can be seen as an IWAE bound with extended latent variable (zki,xmki), whose prior is pθ(xmki|zki)p(zki). The related importance sampling proposal of the ith term is pθ(x m ki|zki)qγ(zki|xoi ), and the observation model is pφ(si|xoi ,xmki)pθ(xoi |zki).
Since all n terms of the sum are IWAE bounds, Theorem 1 from Burda et al. (2016) directly gives the monotonicity property:
L1(θ, φ, γ) ≤ . . . ≤ LK(θ, φ, γ). (23)
Regarding convergence of the bound to the true likelihood, we can use Theorem 3 of Domke & Sheldon (2018) for each term of the sum to get the following result. Theorem. Assuming that, for all i ∈ {1, ..., n},
• there exists αi > 0 such that E [ |w1i − pθ,φ(xoi , si)|2+αi ] <∞,
• lim supK−→∞ E [ K/(w1i + ...+ wKi) ] <∞,
the not-MIWAE bound converges to the true likelihood at rate 1/K:
`(θ, φ)− LK(θ, φ, γ) ∼ K→∞
1
K n∑ i=1 Var[w1i] 2pθ,φ(x o i , si) 2 . (24)
E VARYING MISSING RATE (UCI)
The UCI experiments use a self-masking missing process in half the features: when the feature value is higher than the feature mean it is set to missing. In order to investigate varying missing rates we change the cutoff point from the mean to the mean plus an offset. The offsets used are {0, 0.25, 0.5, 0.75, 1.0}, so that the largest cutoff point will be the mean plus one standard deviation. Increasing the cutoff point further results in mainly imputing outliers. Results for PPCA and notMIWAE PPCA using the agnostic missing model are seen in figure 5 and using the self-masking model with known sign of the weights are seen in figure 6. Figure 7 shows the results for MIWAE and not-MIWAE using self-masking with known sign of the weights.

Review summaries:
This paper proposes an approach to training deep latent variable models on data that is missing not at random. To learn the parameters of deep latent variable models, the paper adopts importance-weighted variational inference techniques. Experiments on a variety of datasets show that the proposed approach is effective by explicitly modeling missing not at random data.