{'fraction_non_alphanumeric': 0.03853705310259806, 'fraction_numerical': 0.01315008447102872, 'mean_word_length': 4.593538500830035, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unsupervised Imputation of Non-ignorably Missing Data Using Importance-Weighted Autoencoders

Deep Learning (DL) methods have dramatically increased in popularity in recent years. While its initial success was demonstrated in the classification and manipulation of image data, there has been significant growth in the application of DL methods to problems in the biomedical sciences. However, the greater prevalence and complexity of missing data in biomedical datasets present significant challenges for DL methods. Here, we provide a formal treatment of missing data in the context of Variational Autoencoders (VAEs), a popular unsupervised DL architecture commonly utilized for dimension reduction, imputation, and learning latent representations of complex data. We propose a new VAE architecture, NIMIWAE, that is one of the first to flexibly account for both ignorable and non-ignorable patterns of missingness in input features at training time. Following training, samples can be drawn from the approximate posterior distribution of the missing data can be used for multiple imputation, facilitating downstream analyses on high dimensional incomplete datasets. We demonstrate through statistical simulation that our method outperforms existing approaches for unsupervised learning tasks and imputation accuracy. We conclude with a case study of an EHR dataset pertaining to 12,000 ICU patients containing a large number of diagnostic measurements and clinical outcomes, where many features are only partially observed.


Introduction
Deep Learning (DL) methods have dramatically increased in popularity in recent years.
While its initial success was demonstrated in the classification and manipulation of image data, there has been significant growth in the application of DL methods to an array of problems in the biomedical sciences (Shickel et al., 2018;Lopez et al., 2018). The presence of complex interactions among high-dimensional features in modern biomedical data has motivated the use of Variational Autoencoders (VAEs), a DL architecture commonly utilized for unsupervised learning tasks such as dimension reduction, representational learning, and generation of synthetic data mimicking real input data, which may be unavailable due to patient confidentiality (Shickel et al., 2018). In prior evaluations, VAEs have shown tremendous performance in data generation and representation learning (Kingma and Welling, 2019).
Some of the useful characteristics of the VAE can be derived from its underlying architecture. This structure is similar to that of the Autoencoder (AE), which consists of two major parts. The first part is the encoder neural network, which in the AE deterministically transforms the input data into a lower-dimensional set of factors that is assumed to capture the salient qualities of the data. The second is the decoder, which is a separate neural network that attempts to reconstruct the original data from the lower dimensional space itself (Tschannen et al., 2018). The quality of the dimension reduction and reconstruction is measured by the reconstruction error, which is measured by the difference in the input data and its reconstructed version. This lower dimensional space may be used to define structure in the underlying set of input data (Wang et al., 2016), similar to PCA. However, the encoder and decoder are often represented by feed forward neural networks, allowing for complex nonlinear interactions between features to be captured in both (LeCun et al., 2015). The VAE additionally imposes a probabilistic assumption on the input data and the lower-dimensional (or "latent") space (Doersch, 2016), where explicit encoder and decoder distributions are defined, and neural networks are used to output the parameters of each distribution (Kingma and Welling, 2019). The learned decoder allows the VAE to also generate synthetic data that has high fidelity to the original training data, given the learned latent space.
However, accounting for the prevalence and complexity of missing data in modern biomedical datasets presents a significant challenge for training VAEs and other DL architectures.
Existing DL methods that claim to handle missing data either require pre-training on studies with fully observed data, utilize heuristics to replace missing observations, or cannot handle more complicated forms of missingness, such as Missing at Random (MAR), or Missing Not a Random (MNAR) (Li et al., 2020;Strauss and Oliva, 2021). Intuitively, when observations with missing values differ systematically from those without missingness, results from standard approaches that do not properly account for this fact may no longer be accurate. For example, the common presence of missing data across patient records in electronic health record data has been shown to present a significant barrier to the generalizability and applicability of deep learning methods (Wells et al., 2013). EHR data also typically contain many features pertaining to a large number of patients (Ross et al., 2014), and such features may show complex interactions with each other as well as with clinical outcomes of interest (Vinzamuri and Reddy, 2013), compounding the impact of the missingness of these features. For these reasons, the ability to flexibly account for different patterns of missing data in modern biomedical datasets, especially in the more difficult nonignorable missingness setting, would have great utility.
In this manuscript, we introduce a novel deep learning architecture called Non-Ignorably Missing Importance-weighted Autoencoder (NIMIWAE) that treats missing observations as latent variables in the VAE framework using Importance-Weighted Autoencoders (IWAEs).
Here we model the probability of missingness using a feed forward neural network, facilitating more flexible handling of MNAR patterns of missingness across a large set of available features. Through simulations, we show that proper modelling of the missingness mechanism increases the accuracy of missing data imputation, as well as downstream coefficient estimation using these imputed datasets. Lastly, we show that in the Physionet 2012 Challenge EHR dataset, accounting for MNAR missingness results in significant differences in downstream fitted models to predict in-hospital mortality. More generally, NIMIWAE is a flexible framework within the VAE family to handle multiple patterns of missingness commonly found in complex biomedical datasets.

Variational Autoencoder
Let X be an n × p data matrix, where x i denotes the observation vector pertaining to the i th observation, i = 1, . . . , n, and x ij denotes the value of the j th feature in this vector, j = 1, . . . , p. In a VAE, we assume x 1 , . . . , x n are i.i.d. samples from a multivariate p.d.f or "generative model" p ψ (X|Z). Here, Z is an n × d matrix, such that Z = {z 1 , · · · , z n } and z i is a latent vector of length d pertaining to the i th sample (Kingma and Welling, 2019).
Typically it is assumed that d p, such that Z constitutes a lower-dimensional representation of the original data X. The parameters ψ and conditional distribution p ψ (X|Z) indicate how the observed data X may be generated from Z. In this manner, a VAE aims to learn accurate representations of high-dimensional data, and may be used to generate synthetic data with similar qualities as its training data. These aspects are also aided through the use of embedded deep learning neural networks, for example within p ψ (X|Z), which also facilitate its applicability to larger dimensions and complex datasets.

Objective Function
Since ψ is unknown, learning is performed by maximizing the marginal log-likelihood of X with respect to ψ, where we denote this marginal log-likelihood as log p ψ (X) = log p ψ (X, Z)dZ = log p ψ (X|Z)p(Z)dZ. However, due to the integral involved, this quantity is often intractable and is difficult to maximize directly. Therefore, VAE's alternatively optimize the so-called "Evidence Lower Bound" (ELBO), which has the following form (Kingma and Welling, 2019): Here, L ELBO (θ, ψ) denotes the ELBO such that L ELBO (θ, ψ) log p ψ (X). Also letL ELBO K (θ, ψ) denote the empirical approximation to Eq. (1) computed by Monte Carlo integration, such that L ELBO (θ, ψ) ≈L ELBO K (θ, ψ) andZ 1 , . . . ,Z K are K samples drawn from q θ (Z|X), the variational approximation of the true but intractable posterior p ψ (Z|X), also called the "recognition model". Furthermore, denote f ψ (Z) and g θ (X) as the decoder and encoder feed forward neural networks of the VAE, where ψ and θ are the sets of weights and biases pertaining to each of these neural networks, respectively. Given Z, f ψ (Z) outputs the distributional parameters pertaining to p ψ (X|Z). Given X, g θ (X) outputs the distributional parameters for q θ (Z|X).
In variational inference, q θ (Z|X) is constrained to be from a class of simple distributions, or "variational family", to obtain the best candidate from within that class to approximate p ψ (Z|X). Variational inference is usually used in tandem with amortization of the parameters where the neural network parameters are shared across observations (Gershman and Goodman, 2014), allowing for stochastic gradient descent (SGD) to be used for optimization of Eq.
(2) (Kingma and Welling, 2019). In practice, both q θ (Z|X) and p(Z) are typically assumed to have simple forms, such as multivariate Gaussians with diagonal covariance structures, and q θ (Z|X) is commonly assumed to be factorizable, such that q θ (Z|X) = n i=1 q θ (z i |x i ) (Kingma and Welling, 2019).

Estimation Procedure and Use Cases
Let (θ (t) ,ψ (t) ) be the estimates of (θ, ψ) at update (or iteration) t. For t = 0, these values are often initialized to small values centered around 0, although other initialization schemes may be used Murphy, 2016). Each subsequent update t 1 consists of two general steps to maximize L(θ, ψ). First, K samples are drawn from qθ (t) (Z|X) to compute the quantity in Eq. (2), conditional on (θ (t) ,ψ (t) ). Then, the so-called "reparametrization trick" is utilized to facilitate the calculation of gradients of this approximation to obtain (θ t+1 ,ψ t+1 ) using stochastic gradient descent (Kingma and Welling, 2019). This procedure may be repeated for a fixed number of iterations, or may be terminated early via pre-specified early stop criteria . Kingma and Welling (2019) provides additional details on the maximization procedure for VAEs. The networks f ψ (Z) and g θ (X) also allow the VAE to capture complex and non-linear relationships between features when outputting the distributional parameters for the generative and recognition models, respectively, through the inclusion of hidden layers in each network. The number of hidden layers and nodes per layer for each network are commonly determined via hyperparameter tuning.
After model fitting, the VAE has several useful features. First, synthetic data can be generated by sampling from the learned generative model pψ(X|Z) after drawing a sample from qθ(Z|X) . Second, the posterior modes in the latent space Z can be determined from qθ(Z|X). These posterior modes may be used for purposes such as clustering or substructure discovery (Lim et al., 2020). In some applications, pψ(X|Z) may also be used to directly perform imputation , however the statistical properties of this procedure have not been thoroughly discussed in prior work.

Importance-Weighted Autoencoder
The IWAE (Burda et al., 2015) is a generalization of the standard VAE, where the resulting IWAE bound, corresponding to the ELBO in Eq. (1), can be written as . Then, to compute its empirical approximationL IW AE K , typically only one sample is drawn for each Z k in Eq. (4). For K > 1, Burda et al. (2015) showed that log p(X) L IW AE Thus, the IWAE bound more closely approximates the true marginal log likelihood when K > 1 (Cremer et al., 2017), at the cost of greater computational burden. If K = 1, L IW AE 1 = L V AE , and the IWAE corresponds exactly to the standard VAE. In this way, the IWAE can be considered to be part of the VAE family, and we refer to methods that use either VAEs or IWAEs broadly as "VAE methods". A visualization of the workflow for an IWAE (without missing data) can be found in Web Appendix A.
VAE methods have shown excellent performance in representation learning on many types of data. However, the presence of missingness in X presents significant challenges to the above modeling procedures and the application of VAE methods in general.

Missing Data
In this section, we first introduce notation for missing data and review the different mechanisms of missingness, as described in the statistical literature. Let the data be factored such that X = {X o , X m }, with X o denoting the observed values and X m denoting the missing values. For each observation vector x i , denote x o i and x m i respectively to be the observed and missing features of x i . Also, let R be a matrix of the same dimension as X, with entries r ij = I(x ij is observed) for the i th observation and j th feature, where I(·) denotes the indicator function. In this way, R is the "mask" matrix pertaining to X, such that x o i = {x ij : r ij = 1} and x m i = {x ij : r ij = 0} for all i = 1, . . . , n and j = 1, . . . , p.
Missingness was classified into three major categories, or mechanisms, in the seminal work by Little and Rubin (2002). These mechanisms are missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR), and they satisfy the fol- Here, φ denotes the unknown parameters pertaining to the missingess model p(r i |x i , z i , φ), where r i = {r i1 , . . . , r ip }. We discuss forms of this model in Section 2.3.2. In the presence of missingness mask R, the marginal loglikelihood may be written as We may factor p ψ,φ (X o , X m , Z, R) using the selection model factorization (Diggle and Kenward, 1994), which is written as

Ignorable Missingness
In a likelihood-based analysis under either MCAR or MAR, the missingness mechanism is considered to be "ignorable" such that the missingness mechanism need not be explicitly modelled in these cases (Rubin, 1976;Little and Rubin, 2002). Under ignorable missingness, the left hand side of Eq. (5) can be separated into log need not be specified because inference on the parameters of interest pertaining to p ψ (X o ) is independent of p φ (R|X o ). Then, one aims to maximize the quantity This quantity can be bounded below exactly as in Section 2.1.1, conditioning on just the observed data X o , rather than the full data X. Existing methods typically take advantage of this simplification, and are shown to perform well under ignorable missingness. Details for these methods are given in Web Appendix A.

Non-ignorable Missingness
In contrast, non-ignorable (or MNAR) missingness refers to the case where the missingness can be dependent on any unobserved values, including the missing entries x m i . MNAR missingness can also be dependent on x o i as well as latent values like Z, and thus MNAR represents the most general and difficult case of missingness in practice. Here, we assume that R is independent of Z, as conditioning on such latent factors may be computationally redundant based on the assumed data generating process (Ibrahim, 2001). In this setting, the missingness typically requires specification of a model for the missingness p(R|X, φ) (Stubbendick and Ibrahim, 2003). Current VAE methods are only able to handle MCAR or MAR missingness, and there is no method to properly deal with the more difficult MNAR case. This issue is especially problematic because missingness in many real world applications have been posited to be non-ignorable (Beaulieu-Jones and and, 2016; O'Shea, 2019).
There have been a number of ways to specify p(R|X, φ) in statistical literature. For example, Diggle and Kenward (1994) proposes a binomial model for the missing data mechanism: where j m = 1, . . . , p miss indexes the p miss features in X that contain missingness, φ jm is the sets of coefficients pertaining to j th m missingness model, and p(r ijm = 1|x i , φ jm ) can be modeled straightforwardly by a logistic regression model, such that where φ 0jm is the intercept of the j th m missingness model, φ 1jm = {φ 1,jm,1 , . . . , φ 1,jm,p obs } T is a p obs × 1 vector of coefficients pertaining to the fully-observed features, and φ 2jm = {φ 2,jm,1 , . . . , φ 2,jm,p miss } T is a p miss ×1 vector of coefficients pertaining to the missing features.

NIMIWAE: IWAE with Nonignorable Missingness
We now propose a novel method to perform statistical learning and imputation using an IWAE in the presence of missing data (NIMIWAE), assuming missingness is nonignorable.
We later show how this model can be simplified when missingness is assumed to be ignorable (IMIWAE). First, we specify a general form of the lower bound, in which we utilize the general IWAE framework to form a tighter bound on the marginal log-likelihood than the VAE ELBO. Let us define q θ (Z, X m ) as the variational joint posterior pertaining to (Z, X m ).
We can factor this variational joint posterior as q θ (Z, pertaining to the missing features, where each X m k has dimensionality p miss . Similar to traditional VAEs, we utilize the class of factorized variational posteriors, except now q θ (Z, Then, denoting z ik and x m ik as the i th observation vectors of Z k and X m k , respectively, we have z i1 , . . . , z iK The form of the lower bound, which we call the NonIgnorably Missing Importance-Weighted Auto Encoder bound, or "NIMIWAE bound", is derived as follows: As explained in Section 2.3, we use the selection model factorization of the joint distribution Here, ψ denotes the weights and biases of the encoder and decoder neural networks, and φ denotes the weights and biases of the missingness network that learns the parameters of the missingness model. Applying the above factorizations to Eq. (8), and estimating the expectations in Eq. (8) by sampling from q θ (z i , x m i ) we obtain the estimate of the NIMIWAE bound: where {z ik ,x m ik } are samples of {z i , x m i } that are drawn via ancestral sampling (Bishop, 2006) , respectively, and the NIMIWAE bound is optimized using the Adam optimizer (Kingma and Ba, 2014). In NIMIWAE, we have four neural . The quantity h φ (x i ), which we call the "missingness network", outputs the distributional parameters to the missingness model p φ (r i |x o i ,x m ik ) specifically given in Eq. (9). Furthermore, by omitting this network and its contribution to the NIMIWAE bound altogether, one can attain an ignorably-missing version of the NIMIWAE method (IMIWAE), which would be more suitable for use when assuming MCAR and MAR missingness. We explore the empirical performance of each of these models under misspecification of the missingness mechanism in Section 3. An illustration of h φ (x i ) is given in Figure 1.
Historically, the set of features for the p miss logistic regression models from Eq. (7) need to be carefully pre-specified, usually based upon prior information (Little and Rubin, 2002).
Prior work has shown that overparameterization of the missingness model can lead to identifiability issues and divergence in EM-based maximization procedures (Ibrahim and Molenberghs, 2009). Our proposed method allows users to similarly pre-specify a subset of features in the missingness network. When such information is not available, it is unclear how such overparameterization will affect model fitting and performance in this setting. Therefore, in Section 3, we evaluate empirically the use of all p features in the missingness network when assuming MNAR missingness using statistical simulation.
Additional details regarding the NIMIWAE algorithm are outlined in Web Appendix A.
[ We first note that Then, we may estimate this integral by self-normalized importance sampling. We utilize , and consider p ψ,φ (x o i , r i ) as some unknown constant. Then, we define the quantities . . , K, with 1 sample drawn from the variational posterior of each latent variable z ik and x m ik to compute s ik , where w ik is defined as standardized "importance weights" . Using these weights we may estimate Then, the process can be repeated for each observation i = 1, . . . , n.
In the MCAR or MAR case, one can similarly estimate E[x m i |x o i ] using the fitted IMIWAE model. By following a similar derivation using the proposal density , with w ik defined as before, but with a slightly different form for s ik : .
Given these weights w ik , we may now also produce Q multiply-imputed datasets using the sampling importance resampling (SIR) algorithm (Smith and Gelfand, 1992). We first construct a set of K candidate draws and corresponding weights using the procedure described above (default K = 10 × Q draws), and then perform a weighted resample of size Q with replacement from this set of draws to obtain approximate draws from p ψ (x m i |x o i , r i ), for each observation i = 1, . . . , n. We may then use techniques such as "Rubin's rules" (Rubin, 2004) to pool the estimates obtained from a candidate regression model fit on each imputed dataset, and obtain pooled point estimates and standard errors that account for the uncertainty due to the imputation. In our analyses, we used NIMIWAE to construct Q = 50 multiply-imputed datasets by drawing K = 500 times from q θ 2 (x m i |z i , x o i , r i ) for each i = 1, . . . , n after the model was trained.

Simulated Data
We utilize statistical simulation to evaluate the imputation performance of our proposed NIMIWAE and IMIWAE methods and under the assumption of MCAR, MAR, and MNAR missingness. We note that for each these simulations, we utilize all p features in NIMIWAE's missingness network. We also compared this performance to state-of-the-art missing data methods in machine learning that claim to handle ignorable missingness patterns: HIVAE , VAEAC , MIWAE , in addition to the popular MICE method (Van Buuren and Groothuis-Oudshoorn, 2011) and a naïve mean imputation method. We also included MissForest (Stekhoven and Buhlmann, 2011) in analyses where the model could be fit in a CPU with 32 GB of memory. For all simulations, we divided the full data into training and validation sets with ratio 8:2. For methods that require hyperparameter tuning, each method was trained on the training set for a given set of hyperparameters, the performance of this model was then evaluated on the validation set, and finally the training set was imputed using the model pertaining to the optimal set of hyperparameters (based up on validation set performnance). For methods that required no hyperparameter tuning, we directly imputed just the training set, for consistency across all methods. In Sections 3.1.1-3.1.2, we evaluate performance on fully synthetic data.
Then, in Section 3.1.3, we simulate missingness into existing datasets from the UCI machine learning repository to preserve non-linearity and interactions between features previously observed. The simulation setup and performance criteria are described in the subsequent sections.

Simulation Setup
We first evaluate the performance of each method on completely synthetic data. Here we assume X is generated such that X = ZW + B, where W and B and are matrices of dimensions d × p and n × p, respectively, Z ∼ N d (0, I), W lj ∼ N (0, 0.5), and B ij ∼ N (0, 1) for i = 1, . . . , n, j = 1, . . . , p, and l = 1, . . . , d.
We then simulate the missingness mask matrix R such that 30% of features are partially observed, and 50% of the observations for each of these features are missing. We generate r ij from the Bernoulli distribution with probability equal to p(r ijm = 1|x i , φ), such that Here, we assume that j m = 1, . . . , p miss index the missing features, φ 1 = {φ 11 , . . . , φ 1,p obs } are the coefficients pertaining to the fully observed features, and φ 2 = {φ 21 , . . . , φ 2,p miss } are those pertaining to the partially observed features, and p obs and p miss are the total number of features that are fully and partially observed, respectively. We set p miss = 0.3 * p and p obs = p − p miss . We drew nonzero values of φ 1 and φ 2 from the log-normal distribution with mean µ φ = 5, with log standard deviation To evaluate the impact of the misspecification of the missingness mechanism on model performance, r ijm was simulated under each mechanism as follows: (1) MCAR: {φ 1 , φ 2 } = 0, (2) MAR: Same as MCAR except φ 1jo = 0 for one randomly selected completely-observed feature j o where j o = p miss + 1, . . . , p, and (3) MNAR: Same as MCAR except φ 2jm = 0. In this way, for each MAR or MNAR feature, the missingness is dependent on just one feature.
In each case, we used φ 0 to adjust the overall expected rate of missingness of each feature to approximately 50%.
Lastly, we simulated a binary response variable assuming Pr(y = 1|X) = Sigmoid(β 0 + Xβ), where y is a binary response variable, β = {β 1 , . . . , β p } are the set of regression coefficients, and β 0 is the intercept. In many applications, it is of interest to use the features in X to predict some outcome variable y when X is only partially observed. Therefore, the ability accurately estimate β is also of importance in the presence of missingness.
Multiple imputation methods like MICE can be used to perform coefficient estimation by using Rubin's rules to pool the coefficient estimates from logistic regression models fitted on each individual multiply-imputed dataset. We similarly perform coefficient estimation using multiply-imputed datasets from the SIR algorithm within NIMIWAE, allowing for direct comparisons in coefficient estimation with MICE in estimating β.
We vary n, p, and d such that n = {10, 000, 100, 000}, p = {25, 100} features, and d = {2, 8}. We simulated 5 datasets per simulation condition, spanning various missingness mechanisms and values for {n, p, d}. We fix the values of β at 0.25 for each feature, and adjusted β 0 to ensure equal proportions for the two binary classes in Y. We measured imputation performance by calculating the average L1 distance between true and imputed masked values in X. LettingX m denote the imputed masked values of the true X m values of the missing entries, we denote the average L1 distance is simply |X m −X m | N miss , where N miss is the total number of missing entries in the dataset. To assess the ability of multiple imputation methods in performing coefficient estimation, we reported the percent bias (PB) of these pooled estimates compared to the truth, averaged across the p features, i.e. P B = 1 [ Figure 2 about here.] Figure 2 shows the results pertaining to n = 100, 000. [ Figure 3 about here.] Figure 3 shows the results pertaining to n = 10, 000 utilizing the alternative initialization method described in Section 2.4, due to the smaller sample size in this setting. We see that NIMIWAE performs best in imputing MNAR values, and in estimating the coefficients under MNAR missingness, and still performs comparably well to other methods in the MCAR and MAR cases. The results of the analyses on n = 10, 000 simulations with default initialization can be found in Web Appendix B. We found that the default method initialized weights for missing features too small for the network to recover in the MNAR case, under the lower sample size and high dimensionality setting. Based upon these results, we recommend the alternative weight initialization for smaller sample sizes (n 10, 000)

Simulation Results
with large dimensionality (p 100), while the default initialization may suffice when the sample size is large, or if the dimensionality of the data is smaller. Alternatively, one may greatly improve imputation performance by narrowing down the features that are input into the missingness network using some prior knowledge or assumptions on the mechanism of missingness. In practical applications, it is unknown which mechanism may truly underly the data, and, as usual, sensitivity analyses, where one varies the assumed mechanism, is still important. For example, imputing data via IMIWAE may be helpful and be more efficient under MCAR or MAR than via NIMIWAE.
Overall, these simulations confirm that existing methods can impute MCAR and MAR missingness with a reasonable degree of accuracy, but they break down under MNAR missingness. Our NIMIWAE method is able to adequately impute MNAR missing values, while still imputing MCAR and MAR missing values with a comparable degree of accuracy as existing methods geared specifically towards the ignorable mechanisms of missingness. Additionally, we note that MissForest was able to be run to the n = 10, 000 and p = 25 simulation case only due to memory constraints, and yielded very poor imputation performance in the MNAR case, like the other ignorably-missing methods. We also show empirically that using all p features in NIMIWAE's missingness network can yield reasonable performance, especially when the number of samples is large. We postulate that this may be due to the fact that neural networks have been shown to generalize well despite severe overparameterization (Poggio et al., 2020). Still, the specification of a smaller model that is closer to the truth may improve the accuracy of imputations, especially under smaller sample sizes (Du et al., 2021).

UCI Machine Learning Datasets
Next, we analyzed how accurately these methods can impute missing values that may be present in real data. Since we are unable to obtain the true missing values in real datasets, we selected 6 different completely-observed datasets from the UCI Machine Learning Repository, and simulated missingness according to each mechanism as in the fully synthetic datasets.
Here, we masked half of the features, such that p miss = f loor(0.5 * p). Additional details on these datasets and how they may be obtained can be found in Web Appendix B, and the values of the hyperparameters that were tuned can be found in Web Appendix C. [

Physionet 2012 Challenge Dataset
Finally, we analyzed the Physionet 2012 Challenge data using each of the compared methods. This data consisted of 114 features of 12,000 patients after pre-processing, and the details of the acquisition and preparation of the data can be found in Web Appendix C. We performed a qualitative analysis of each imputed dataset, highlighting differences between the results of downstream regression models fitted on mortality, assuming non-ignorable (NIMIWAE) and ignorable missingness (Ibrahim et al., 2005;Ibrahim and Molenberghs, 2009). This is because, in contrast to our simulations, the true values of the missing entries are not available to directly assess imputation performance. Additionally, the missingness mechanism itself is generally not "testable" by the observed data in practice (Ibrahim et al., 1999). Here, we used the alternative initialization scheme of NIMIWAE, due to the limited sample size and large dimension of the data. Additionally, we added "supervised" versions of NIMIWAE, IMIWAE, and MICE in our comparisons, where the response of interest (mortality) is included with the features of the data during training and imputation, such that the response variable may inform the multiple imputations. After model fitting, we performed Rubin's rules to pool the estimates of the fitted models on the multiply-imputed datasets, and compute standard errors across imputations.
[ Figure 4 about here.] Based on the imputed datasets from these methods, we fit a logistic regression model with post-baseline mortality as the binary response, with the p = 114 baseline features of this dataset as the covariates. We report details of the covariates with the top 10 largest effects on mortality when using the multiply imputed datasets from the supervised non-ignorably missing NIMIWAE model in Table 4. We found some significant differences in the results of the logistic regression based on the different imputed datasets. For example, we found that N IM IW AE sup uncovered a stronger effect of the variables T emp_highest, K_last and Gender compared to other methods. These factors have been studied for their association with mortality in ICU patients. Specifically, gender has been studied for having a potentially significant effect on mortality in critically-ill patients (Mahmood et al., 2012;Larsson et al., 2019). Additionally, body temperature (Schell-Chaple et al., 2015) and irregular potassium levels (Tongyoo et al., 2018) have both been known to be associated with an increased risk of in-hospital mortality.
[ Figure 5 about here.] values of FiO2 that were significantly higher than the observed mean, and values of RespRate that were significantly lower than the observed mean. Some studies have shown that abnormally high or low respiratory rates may be associated with higher mortality in critically-ill patients (Strauß et al., 2014;Ljunggren et al., 2016), suggesting that the missing values of RespRate_last, which were more prevalent in deceased patients, may have truly been further away from the normal range of 12-20. Additionally, Esteban (2002) found that higher levels of administered FiO2 were associated linearly with higher mortality. Thus, the missing values of F iO2_last, which were more prevalent in surviving patients, may have been lower than the observed values. This suggests that the mechanism of missingness of F iO2_last and RespRate_last may be non-ignorable, since N IM IW AE imputed more realistic values for these variables according to the mortality rates of patients with missingness in these variables.

Discussion
In this paper we introduce NIMIWAE, one of the first methods to handle up to MNAR patterns of missingness in the VAE/IWAE class of methods, to address complex patterns of missingess observed in the Physionet EHR data. Using statistical simulations, we show that NIMIWAE performs well in imputing missing features under MNAR, and has reasonable performance under the MCAR and MAR settings. Performance in imputing MCAR and MAR missingness can be further improved in NIMIWAE by using the ignorable version of this model (IMIWAE), where we omit the missingness network. We also found that the results of the analysis on the Physionet data are highly dependent on the choice of missingness model, which specifies the assumption of the missingness mechanism. However, NIMIWAE is able to impute missing values well in simulations regardless of the underlying missingness mechanism, flexibly modelling the mechanism using a deeply-learned neural network. The NIMIWAE-imputed dataset resulted in more realistic imputed values with respect to what we may expect in the Physionet Challenge patients, since NIMIWAE takes into account possible non-ignorable missingness in the data. Additionally, the NIMIWAE architecture learns a lower-dimensional representation of the data, which can be used for tasks such as patient subgroup identification or visualization of data.
Learning algorithms that can be applied to EHR data like the Physionet Challenge dataset can be valuable tools that clinicians can use to aid decisions in hospital settings and understand patterns within these health records. For example, properly handling missingness in EHRs when imputing the missing entries can improve the performance of prediction algorithms that can assess risk of death or other outcomes of interest, like disease. Informative missingness is a common problem in analyzing EHR data, and accounting for such missingness can be helpful in obtaining accurate, unbiased estimates of the true missing values. We note that although we have used our NIMIWAE method primarily to analyze the Physionet 2012 Challenge dataset, it can more generally be applied to settings where one wishes to train a VAE when missingness is present among input features.       (1994) (as given in Section 2.3.2 of the main text) when h φ (x i ) contains no hidden layers, and the selection of additional hidden layers can allow for the capturing of more complex effects.

Medical
Under simple distributional assumptions of q θ 2 (x m i |z i , r i , x o i ), the sampling step in Step (4) is similar to the sampling of the latent variable z i in Step (2), and both can be accomplished using the reparametrization trick (Kingma & Welling 2013).

NIMIWAE: Initialization, Early Stop, and Hyperparameter Tuning
Initialization of weights and biases in deep learning architectures can significantly affect the trained model, especially in datasets with smaller sample sizes. By default, NIMIWAE uses the semi-orthogonal matrix initialization  for ψ, θ, and φ. Alternatively, we initialize weights pertaining to the missing features in the input layer of the missingness network using values drawn from a Uniform(-2,2) distribution. This is done in order to draw larger initial values of the effects of missing variables on the missingness, in order to encourage the network to learn nonzero effects of these missing features. Empirically, we found that this alternative initialization helps the network impute more accurate values in smaller sample size settings, particularly under MNAR missingness, while maintaining similar performance in the MCAR and MAR cases.
We also incorporate an early stop criterion In this way, training continues as long as the improvement in the estimated lower bound in the validation set is greater than εL opt , but if not, training is allowed to run for a set number of epochs before early stopping. This leeway (called "patience") is allowed due to the properties of SGD, which may cause L (τ ) to fluctuate, especially for smaller τ . If E (τ ) = patience, we save the optimal model pertaining to L opt , and terminate training. Here, we set patience = 50 and ε = 0.0001.
In order to find the optimal combination of values for these hyperparameters, we perform an extensive grid search of all combinations of prespecified values for these quantities. One may also choose to tune the number of latent variables K. In this paper, we fix K = 5 during training for computational ease. 6

Web Appendix B: Additional Simulations and Details
Additional Simulations Web Figure 2 shows results from the simulation setup with n = 10, 000, and Web Figure   3 shows results from the alternate setup with n = 100, 000, varying p = {25, 100} and d = {2, 8} in each case.
As in all deep learning algorithms, NIMIWAE performs best under larger sample sizes, and this aspect is especially important given the additional task in NIMIWAE of narrowing down relevant features in a potentially overspecified missingness model. When a large number of samples is not readily available, and the dimensionality of the data is very large, one can either utilize an alternative method of weight initialization for the missing features in the missingness network, as seen in Section 3.1 of the main text, or narrow down the number of features that are included in NIMIWAE's missingness network using some prior knowledge on the dependencies of the features with the missingness. We found that under the large sample size setting (n = 100, 000), imputation performance was generally robust to the method of initialization, as shown in Web Figure 3. We additionally found that having fewer unnecessary features in the missingness model greatly increased the performance of both imputation and inference under MNAR, as one would expect. However, omitting a feature that is truly relevant to the missingness may cause very poor performance in both imputation and downstream coefficient estimation, and in practice, determining which features are relevant to the missingness may be nontrivial. One may utilize a priori knowledge about the data in order to inform a decision on which features are significant in the missingness model, but such decisions may not be able to be made with a high levels of certainty.
Web Figure 2: Average L1 distance between true and imputed values for missing entries stratified by d (left), and percent bias of pooled coefficient estimates. NIMIWAE outperforms all methods in imputing missing values that were simulated to be MNAR when p = 25, but performs poorly under p = 100. Here, n = 10, 000, µ φ = 5, and error bars show the variability of each metric across the 5 reps. NIMIWAE struggles to handle the difficult MNAR missingness pattern when the sample size is smaller, while the dimensionality of the data is large.
Web Figure 3: Average L1 distance between true and imputed values for missing entries stratified by d (left), and percent bias of pooled coefficient estimates. NIMIWAE was run using the alternative initialization method, as described in Section 2.3 of the main text. Again, NIMIWAE outperforms all methods in imputing missing values that were simulated to be MNAR. Here, n = 100, 000, µ φ = 5, and error bars show the variability of each metric across the 5 reps.

Tuning Hyperparameters
In this section, we summarize the combinations of hyperparameter values that were searched.
For each method and dataset, we searched over combinations of two different values for each hyperparameter that was tuned. We tuned several hyperparameters for each method: number of nodes per hidden layer (h), number of nodes per hidden layer in the missingness network (h r ) , number of hidden layers (nhl), number of hidden layers in the missingness network (nhl r ), dimensionality of latent space (dz), and learning rate (lr). For HI-VAE, nhl was not able to be varied. Therefore, we fixed nhl and tuned the dimensionality of the y latent variable (dy) instead. Hyperparameters nhl r and h r were applicable only to NIMIWAE, so they were tuned only for NIMIWAE. Hyperparameters batch size (bs) and

Web Appendix C: Data
Datasets used for analysis in this study are all publicly available online. Details of how to obtain these datasets are given in this section.

UCI Machine Learning Datasets
The UCI datasets can be obtained from the UCI Machine Learning repository at https: //archive.ics.uci.edu/ml/index.php (Dua & Graff 2017) The datasets included in the analyses in this study are given below: 48 hours. Example features include albumin levels (in g/dL), serum sodium (in mEq/L) and white blood cell count (in cells/nL). For each patient, clinical outcomes such as the SAPS-I score (Gall et al. 1984), SOFA score (Ferreira 2001), length of stay, survival, and in-hospital death were also recorded. The challenge called for learning algorithms to accurately predict the outcome of interest (like in-hospital death) based on such clinical data, facilitating clinical decision-making regarding the need for early intervention in ICU patients (Silva et al. 2012).
The original challenge dataset was divided into training, validation, and test datasets of equal sizes, and missingness was highly prevalent in all three datasets. Table 1 shows the missingness percentage of all features in the dataset, and the number of subjects with at least one non-missing entry recorded for each respective feature.
The data was pre-processed following Johnson et al. (2012), and this procedure, as well as a copy of this data can be found here: https://github.com/alistairewj/challenge2012.
Domain knowledge and distributional assumptions were used in order to convert implausible values, and temporal data was processed by using domain knowledge to break down each feature by summary statistics: first, last, lowest, highest, and median values across the time points. The resulting dataset contained 114 features for each of the 12,000 patients.
In this way, a measurement is missing for a patient only if the diagnostic measure was not taken during any of the 48 time points. We also split the patients with the same 8:2 training-