# Deeply-Learned Generalized Linear Models with Missing Data

David K. Lim
Department of Biostatistics, University of North Carolina at Chapel Hill
and
Naim U. Rashid
Department of Biostatistics, University of North Carolina at Chapel Hill
and
Junier B. Oliva
Department of Computer Science, University of North Carolina at Chapel Hill
and
Joseph G. Ibrahim
Department of Biostatistics, University of North Carolina at Chapel Hill The
authors gratefully acknowledge NIH grants U01-CA274298, P50-CA257911,
P50-CA058223, T32-CA106209, 1R01AA02687901A1, and 1OT2OD032581-02-321, and NSF
grants IIS2133595 and DMS2324394 for funding this research.

###### Abstract

Deep Learning (DL) methods have dramatically increased in popularity in recent
years, with significant growth in their application to various supervised
learning problems. However, the greater prevalence and complexity of missing
data in such datasets present significant challenges for DL methods. Here, we
provide a formal treatment of missing data in the context of deeply learned
generalized linear models, a supervised DL architecture for regression and
classification problems. We propose a new architecture, dlglm, that is one of
the first to be able to flexibly account for both ignorable and non-ignorable
patterns of missingness in input features and response at training time. We
demonstrate through statistical simulation that our method outperforms
existing approaches for supervised learning tasks in the presence of missing
not at random (MNAR) missingness. We conclude with a case study of the Bank
Marketing dataset from the UCI Machine Learning Repository, in which we
predict whether clients subscribed to a product based on phone survey data.
Supplementary materials for this article are available online.

Keywords: missing data, supervised learning, deeply learned glm, MNAR

##  1 Introduction

Deep Learning (DL) methods have been increasingly used in an array of
supervised learning problems in various fields, for example, in the biomedical
sciences (Razzak et al., 2017; Lopez et al., 2018). While a number of deep
learning architectures have been proposed for supervised learning, the feed
forward neural network (FFNN) is commonly used in most architectures. In a
FFNN, sequential non-linear transformations are applied to the values of the
input layer. Each value in the subsequent layer of the FFNN is computed by
applying a non-linear (or “activation”) function to the linear transformation
of the values in the previous layer, outputting a complex non-linear
transformation of the input (Svozil et al., 1997). For example, a FFNN
architecture called the deeply-learned GLM (Tran et al., 2019) has been
applied in the context of supervised learning to describe nonlinear
relationships between the covariates and the response. Due to the large number
of parameters, so optimization is often done via stochastic gradient descent
(Guo & Gelfand, 1990).

However, the common presence of missing data in datasets can hinder the
training and generalizability of supervised deep learning methods (Wells et
al., 2013), where missingness can occur both in the input features and the
response variable. Missingness has commonly been categorized into three
mechanisms: Missing Completely At Random (MCAR), Missing At Random (MAR), and
Missing Not At Random (MNAR) (Rubin, 1976). While a number of methods have
been proposed in the statistical literature to address MNAR missingness in the
regression setting (Ibrahim et al., 2005), such methods often cannot take into
account complex relationships between predictors and response and are not
scalable to higher dimensions (Chen et al., 2019), or have been specifically
designed for unsupervised learning tasks (Lim et al., 2021). Supervised deep
learning is one way to address capture complex relationships between
predictors and response in a scalable manner (Kingma & Welling, 2019), however
it is unclear how best to account for more complex forms of missingness, such
as MAR or MNAR missingness, in this setting.

There have been some recent attempts to perform prediction using deep learning
in the presence of missing features (Ipsen et al., 2021), but such methods
typically assume either MCAR or MAR missingness. Commonly used missing data
methods in supervised deep learning applications, such as mean imputation or
complete case analysis, have historically yielded biased results (Ibrahim &
Molenberghs, 2009). Multiple Imputation by Chained Equations (mice) has also
been widely employed to account for missing data in a supervised learning
setting. However, mice is unable to apply a trained imputation model to handle
missingness that may exist at test time (Hoogland et al., 2020). In addition,
multiple imputation-based methods may not be feasible to apply when the
downstream model is computationally intensive, such as in the setting of
training a deep learning neural network, since one must train the model
separately for each imputed dataset. Lastly, existing approaches to handle MAR
or MCAR missingness when training deep learning models for supervised learning
tasks are currently limited, and have not been sufficiently explored in the
literature.

To address these issues, we present dlglm: a deep generalized linear model
(GLM) for probabilistic supervised learning in the presence of missing input
features and/or response across a variety of missingness patterns. Our
proposed method utilizes variational inference to learn approximate posterior
distributions for the missing variables, and replaces missing entries with
samples from these distributions during maximization. In this way, dlglm can
perform supervised learning in the presence of missingness in both the
features and the response of interest. We also incorporate a model for the
missingness, which can take into account MNAR patterns of missingness, even at
training time. Through neural networks, dlglm is able to model complex non-
linear relationships between the input features and the response, and is
scalable to large quantities and dimensionalities of data. Prediction can be
done seamlessly on fully- or partially-observed samples using the trained
model, without requiring separate imputation of the missing values.

##  2 Methods

Here we first discuss the formulation of the generalized linear model (GLM) in
Section 2.1, and then introduce the deeply-learned GLM in Section 2.2. We then
discuss missingness in the context of GLMs in Section 2.3, and lastly propose
a novel deep learning architecture dlglm in Section 2.4 to fit deeply learned
GLMs in the presence of missingness.

###  2.1 Generalized Linear Models (GLMs)

Let $\mathbf{X}$ be the $n\times p$ matrix of covariates (input features) with
observation vectors $\mathbf{x}_{i}$, where each corresponding entry $x_{ij}$
denotes the value of the $i^{th}$ observation of the $j^{th}$ feature for
$i=1,\ldots,n$ and $j=1,\ldots,p$. Also, let
$\mathbf{Y}=\\{y_{1},\ldots,y_{n}\\}$ be the vector of univariate responses
where $y_{i}$ is the response pertaining to the $i^{th}$ observation. We note
that $y_{i}$ may also be assumed to be multivariate; however, we focus
specifically on the case of univariate response to simplify the discussion,
and discuss extensions to the setting of multivariate response in Section 4.
Then, denote $\boldsymbol{\eta}=\mathbf{X}\boldsymbol{\beta}$, where
$\boldsymbol{\beta}$ is a vector of regression coefficients and
$\boldsymbol{\eta}$ is the linear predictor. Also define
$\boldsymbol{\mu}=\\{\mu_{1},\ldots,\mu_{n}\\}$ with
$\mu_{i}=E(y_{i}|\mathbf{x}_{i})$ and link function $g(\cdot)$ such that
$g(\mu_{i})=\eta_{i}=\mathbf{x}_{i}\boldsymbol{\beta}$. We assume that the
conditional distribution $p(y_{i}|\mathbf{x}_{i})$ is a member of the
exponential family of distributions (McCullagh & Nelder, 2019), such that
$p(y_{i}|\mathbf{x}_{i})$ can be written as

$p(y_{i}|\mathbf{x}_{i})=\text{exp}\left[\frac{y_{i}\Theta_{i}-b(\Theta_{i})}{a(\alpha)}+c(y_{i},\alpha)\right],$

with canonical parameter $\Theta_{i}$, dispersion parameter $\alpha$, and some
functions $a(\cdot)$, $b(\cdot)$, and $c(\cdot)$. Here, we further assume
$g(\cdot)$ is a canonical link function such that $g(\mu_{i})=\Theta_{i}$.
With the appropriate specification of the canonical link $g(\cdot)$ and
variance function $V_{\alpha}(\cdot)$, we obtain the formulation of a GLM.

GLMs were first motivated by the limitations of the traditional linear model,
which imposed strict assumptions of linearity between $\boldsymbol{\mu}$ and
$\mathbf{X}$ and of normality of errors with fixed variance. GLMs instead
utilize specific link and variance functions, allowing for model fitting on
types of response data that may violate these assumptions, such as count or
categorical outcomes, without having to rely on heuristic transformations of
the data (Nelder & Wedderburn, 1972). Typically, GLMs are estimated by
utilizing iteratively re-weighted least squares in lower dimensions (Holland &
Welsch, 1977), with extensions to the higher dimensional case via penalized
likelihood (Friedman et al., 2010).

###  2.2 Deeply Learned GLMs

The traditional GLM assumes $g(\mu_{i})$ is a linear function of
$\mathbf{x}_{i}$, i.e. $g(\mu_{i})=\mathbf{x}_{i}\boldsymbol{\beta}$. In many
modern applications, one may wish to model $g(\mu_{i})$ as a non-linear
function of $\mathbf{x}_{i}$ or capture complex interactions between features
to predict response (Qi & Wu, 2003). In such cases, we may generalize the GLM
to a deeply-learned GLM (Tran et al., 2019) with the following expression:
$g(\mu_{i})=\eta_{i}=h_{\pi}(\mathbf{x}_{i})\boldsymbol{\beta}$, where
$h_{\pi}(\cdot)$ denotes the output of a series of non-linear transformations
applied to the input $\mathbf{X}$ by a neural network, with weights and bias
parameters denoted by $\pi$. In addition, $\eta_{i}$ can alternatively be
expressed $\eta_{i}=s_{\pi,\boldsymbol{\beta}}(\mathbf{x}_{i})$, where
$s_{\pi,\boldsymbol{\beta}}(\cdot)$ is a neural network where
$\boldsymbol{\beta}$ denotes the weights and bias associated with the output
(last) layer of $s_{\pi,\boldsymbol{\beta}}(\cdot)$. This formulation allows
for the traditional interpretation of $\boldsymbol{\beta}$ as the coefficients
pertaining to a transformed version of the input covariates. Figure 1 shows an
illustration of this architecture.

Figure 1: Visualization of a sample deeply-learned GLM architecture
$s_{\pi,\boldsymbol{\beta}}(\mathbf{x}_{i})$. Here, $\pi$ denotes the set of
weights and biases pertaining to the portion of the architecture from the
input layer to the second to last layer (hidden layer 2).
$h_{\pi}(\mathbf{x}_{i})$ is a subset of the entire architecture, such that
$s_{\pi,\boldsymbol{\beta}}(\mathbf{x}_{i})=h_{\pi}(\mathbf{x}_{i})\boldsymbol{\beta}$.
Original artwork of a feed-forward neural network (Dormehl, 2019) was modified
to show deeply-learned GLM architecture.

Let $n_{HL}$ denote the number of hidden layers in
$s_{\pi,\boldsymbol{\beta}}(\cdot)$. We note that if $n_{HL}=0$, then
$h_{\pi}(\mathbf{x}_{i})=\mathbf{x}_{i}$ and
$s_{\boldsymbol{\beta}}(\mathbf{x}_{i})=\mathbf{x}_{i}\boldsymbol{\beta}$,
reducing to the traditional GLM. Deeply learned GLMs and other neural networks
are often maximized using stochastic gradient descent (Bottou, 2012). Details
of this algorithm can be found in Appendix A1 of the supplementary materials.

###  2.3 Missingness in GLMs

Many modern datasets often contain complex forms of missingness (Ghorbani &
Zou, 2018). In GLMs, missingness can exist in either X or Y. Therefore, we
specify three cases of missingness in this context: missing covariates with
fully-observed response (Case x), missing response with fully-observed
covariates (Case y), and missing covariates and missing response (Case xy).
Define $\mathbf{R}=\\{\mathbf{R}^{X},\mathbf{R}^{Y}\\}$ as the “missingness
mask”, which denotes the missingness of $\\{\mathbf{X},\mathbf{Y}\\}$, such
that $\mathbf{R}^{X}$ and $\mathbf{R}^{Y}$ have the same dimension as
$\mathbf{X}$ and $\mathbf{Y}$, respectively, and a value of 1 in $\mathbf{R}$
denotes that the corresponding entry in $\\{\mathbf{X},\mathbf{Y}\\}$ is
observed, while a value of 0 denotes that it is unobserved. Additionally, let
$\mathbf{R}=\\{\mathbf{r}_{1},\ldots,\mathbf{r}_{n}\\}$ with
$\mathbf{r}_{i}=\\{\mathbf{r}_{i}^{X},r_{i}^{Y}\\}=\\{r_{i1}^{X},\ldots,r_{ip}^{X},r_{i}^{Y}\\}$,
$\mathbf{R}^{X}=\\{\mathbf{r}^{X}_{1},\ldots,\mathbf{r}^{X}_{n}\\}$, and
$\mathbf{R}^{Y}=\\{r^{Y}_{1},\ldots,r^{Y}_{n}\\}$ with elements
$\mathbf{r}^{X}_{i}$ and $r^{Y}_{i}$ pertaining to the missingness of the
$i^{th}$ observation of $\mathbf{X}$ and $\mathbf{Y}$, respectively. Then,
$\mathbf{X}$ and $\mathbf{Y}$ can be factored into the unobserved and observed
entries $\\{\mathbf{X}^{m},\mathbf{X}^{o}\\}$ and
$\\{\mathbf{Y}^{m},\mathbf{Y}^{o}\\}$, respectively, such that
$\mathbf{X}^{m}=\\{\mathbf{X}:\mathbf{R}^{X}=0\\}$ with
$\mathbf{x}_{i}^{m}=\\{\mathbf{x}_{i}:\mathbf{r}_{i}^{X}=0\\}$,
$\mathbf{X}^{o}=\\{\mathbf{X}:\mathbf{R}^{X}=1\\}$ with
$\mathbf{x}_{i}^{o}=\\{\mathbf{x}_{i}:\mathbf{r}_{i}^{X}=1\\}$, and
$\mathbf{Y}^{m}=\\{\mathbf{Y}:\mathbf{R}^{Y}=0\\}$ and
$\mathbf{Y}^{o}=\\{\mathbf{Y}:\mathbf{R}^{Y}=1\\}$, with
$y_{i}^{m}=\\{y_{i}:r_{i}^{Y}=0\\}$ and $y_{i}^{o}=\\{y_{i}:r_{i}^{Y}=1\\}$.

Missingness was classified into three primary mechanisms in the seminal work
by Little & Rubin (2002): missing completely at random (MCAR), missing at
random (MAR), and missing not at random (MNAR). They satisfy the following
relations:

  * •

MCAR:
$p(\mathbf{r}_{i}|\mathbf{x}_{i},y_{i},\boldsymbol{\phi})=p(\mathbf{r}_{i}|\boldsymbol{\phi})$

  * •

MAR:
$p(\mathbf{r}_{i}|\mathbf{x}_{i},y_{i},\boldsymbol{\phi})=p(\mathbf{r}_{i}|\mathbf{x}_{i}^{o},y_{i}^{o},\boldsymbol{\phi})$

  * •

MNAR:
$p(\mathbf{r}_{i}|\mathbf{x}_{i},y_{i},\boldsymbol{\phi})=p(\mathbf{r}_{i}|\mathbf{x}_{i}^{o},\mathbf{x}_{i}^{m},y_{i}^{o},y_{i}^{m},\boldsymbol{\phi})$.

Here, $\boldsymbol{\phi}$ denotes the collection of parameters for the model
of the missingness mask $\mathbf{r}_{i}$. In the presence of missingness, the
marginal log-likelihood can generally be written as

$\displaystyle\log
p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{X}^{o},\mathbf{Y}^{o},\mathbf{R})$
$\displaystyle=\log\iint
p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{X},\mathbf{Y},\mathbf{R})d\mathbf{X}^{m}d\mathbf{Y}^{m}$
$\displaystyle=\log\iint
p_{\alpha,\boldsymbol{\beta},\pi}(\mathbf{Y}|\mathbf{X})p_{\psi}(\mathbf{X})p_{\phi}(\mathbf{R}|\mathbf{X},\mathbf{Y})d\mathbf{X}^{m}d\mathbf{Y}^{m},$
(1)

where $\psi$ is a set of parameters associated with the covariate distribution
$p_{\psi}(\mathbf{X})$. We factor
$p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{X},\mathbf{Y},\mathbf{R})$
using the selection model factorization (Diggle & Kenward, 1994).

Under MNAR, it is not possible to remove
$p_{\phi}(\mathbf{R}|\mathbf{X},\mathbf{Y})$ from the integral, since
$\mathbf{R}$ can depend on $\\{\mathbf{X}^{m},\mathbf{Y}^{m}\\}$. Therefore,
MNAR missingness is said to be non-ignorable, because it requires
specification of the so-called “missingness model”
$p(\mathbf{r}_{i}|\mathbf{x}_{i},y_{i},\boldsymbol{\phi})$ (Stubbendick &
Ibrahim, 2003). There are a number of ways to specify this model. For example,
Diggle & Kenward (1994) proposes a binomial model for the missing data
mechanism, which can be written in this setting as

$p(\mathbf{R}|\mathbf{X},\mathbf{Y},\boldsymbol{\phi}_{j_{m}})=\prod_{i=1}^{n}\prod_{j_{m}=1}^{p_{miss}}\left[p(r_{ij_{m}}=1|\mathbf{x}_{i},y_{i},\boldsymbol{\phi}_{j_{m}})\right]^{r_{ij_{m}}}\left[1-p(r_{ij_{m}}=1|\mathbf{x}_{i},y_{i},\boldsymbol{\phi}_{j_{m}})\right]^{1-r_{ij_{m}}},$

where $j_{m}=1,\ldots,p_{miss}$ indexes the $p_{miss}$ features in
$\\{\mathbf{X},\mathbf{Y}\\}$ that contain missingness. Here
$p_{miss}=p_{miss}^{X}+p_{miss}^{Y}$, where $p_{miss}^{X}$ is the total number
of features containing missingness in $\mathbf{X}$, and $p_{miss}^{Y}$ is 1 if
$\mathbf{Y}$ contains missingness (0 otherwise). Also,
$\boldsymbol{\phi}_{j_{m}}$ is the set of coefficients pertaining to the
missingness model of the $j_{m}^{th}$ missing variable, and
$p(r_{ij_{m}}=1|\mathbf{x}_{i},y_{i},\boldsymbol{\phi}_{j_{m}})$ can be
modeled straightforwardly by a logistic regression model, such that

$\text{logit}[p(r_{ij_{m}}=1|\mathbf{x}_{i},y_{i},\boldsymbol{\phi}_{j_{m}})]=\phi_{0j_{m}}+y_{i}\phi_{1j_{m}}+\mathbf{x}_{i}^{o}\boldsymbol{\phi}_{2j_{m}}+\mathbf{x}_{i}^{m}\boldsymbol{\phi}_{3j_{m}},$

where $\phi_{0j_{m}}$ is the intercept of the $j_{m}^{th}$ missingness model,
$\phi_{1j_{m}}$ is the coefficient pertaining to the response variable
$\mathbf{Y}$, and
$\boldsymbol{\phi}_{2j_{m}}=\\{\phi_{2,j_{m},1},\ldots,\phi_{2,j_{m},p_{obs}^{X}}\\}^{T}$
and
$\boldsymbol{\phi}_{3j_{m}}=\\{\phi_{3,j_{m},1},\ldots,\phi_{3,j_{m},p_{miss}^{X}}\\}^{T}$
are the sets of coefficients of the $j_{m}^{th}$ variable’s missingness model
pertaining to the effects of the observed and missing features on the
missingness, respectively, with $p_{obs}^{X}$ and $p_{miss}^{X}$ denoting the
number of completely-observed and partially observed features in $\mathbf{X}$,
respectively. Note that this model assumes independence of $\mathbf{R}$ across
the $p_{miss}$ missing features, such that the missingness of each variable is
conditionally independent of whether any other variable has been observed,
which may or may not be realistic in some settings (Ibrahim et al., 2005).

When missingness is assumed to be MAR or MCAR, the marginal log-likelihood can
be factored as $\log
p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{X}^{o},\mathbf{Y}^{o},\mathbf{R})=\log
p_{\alpha,\boldsymbol{\beta},\pi,\psi}(\mathbf{X}^{o},\mathbf{Y}^{o})+\log
p_{\phi}(\mathbf{R}|\mathbf{X}^{o},\mathbf{Y^{o}})$. In this case, the
quantity $\log p_{\phi}(\mathbf{R}|\mathbf{X}^{o},\mathbf{Y^{o}})$ need not be
specified, since it is independent from the parameters of interest pertaining
to $p_{\alpha,\boldsymbol{\beta},\pi,\psi}(\mathbf{X}^{o},\mathbf{Y}^{o})$.
Therefore, MAR or MCAR missingness is often referred to as “ignorable”
missingness. Equation (1) can then be expressed as

$\log
p_{\alpha,\boldsymbol{\beta},\pi,\psi}(\mathbf{X}^{o},\mathbf{Y}^{o})=\log\iint
p_{\alpha,\boldsymbol{\beta},\pi}(\mathbf{Y}|\mathbf{X})p_{\psi}(\mathbf{X})d\mathbf{X}^{m}d\mathbf{Y}^{m}.$
(2)

###  2.4 Deeply-learned GLM with Missingness (dlglm)

In this section, we propose an algorithm for training deeply-learned GLMs in
the presence of MCAR, MAR, and MNAR missingness. Before discussing this model,
we first discuss the specification of the so-called covariate distribution
$p_{\psi}(\mathbf{X})$ introduced in Equations 1 and 2, which is critical for
maximizing the marginal log-likelihood in either setting. In Sections
2.4.1-2.4.2, we discuss two different models for $p_{\psi}(\mathbf{X})$, and
then in Section 2.4.3 we propose a novel method to handle missingness using a
deeply-learned GLM architecture with an Importance-Weighted Autoencoder (IWAE)
covariate structure. To simplify the discussion, we narrow the scope of our
discussion to the Case x setting, where only $\mathbf{X}$ contains
missingness, but note that the proposed methodology naturally extends to Case
y and Case xy settings as well.

####  2.4.1 Modeling $p_{\psi}(\mathbf{X})$ with known distribution

Given Eq. 1, we must model X with some assumed covariate distribution
$p_{\psi}(\mathbf{X})$. Care must be taken in specifying this distribution, as
improper specification may reduce the accuracy of estimation of the parameters
of interest $\boldsymbol{\beta}$ (Lipsitz & Ibrahim, 1996). For example, we
may assume $p_{\psi}(\mathbf{X})$ follows some known multivariate distribution
such as the multivariate normal distribution, where $\mathbf{X}\sim
N_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ and
$\psi=\\{\boldsymbol{\mu},\boldsymbol{\Sigma}\\}$. Here, $\psi$ can be
optimized jointly with the rest of the parameters
$\\{\alpha,\boldsymbol{\beta},\pi,\phi\\}$ that are involved in the marginal
log-likelihood. However, this assumption may not be applicable in many
instances such as in the case when $\mathbf{X}$ contains mixed data types,
where both continuous and discrete features may be correlated and a joint
distribution may be difficult to specify in closed form. In certain cases, it
may be beneficial to model $p_{\psi}(\mathbf{X})$ flexibly, such that no
strong prior assumptions need to be made on the form of this distribution. To
address this, a sequence of 1-D conditionals have previously been proposed to
model the covariate distribution (Lipsitz & Ibrahim, 1996), but such a model
may be computationally intractable when the number of covariates is very
large.

Once an explicit form for the covariate distribution is specified, one aims to
maximize the marginal log-likelihood, as introduced in Equation (1) in Section
2.3. However, due to the integral involved, this quantity is often intractable
and is difficult to maximize directly, so a lower bound of the marginal log-
likelihood is often maximized instead. The derivation of this lower bound can
be found in Appendix A2 of the supplementary materials.

####  2.4.2 Modelling $p_{\psi}(\mathbf{X})$ with Variational and Importance-
Weighted Autoencoders

Alternatively, one can approximately learn $p_{\psi}(\mathbf{X})$ from the
training data by using an IWAE neural network architecture. In this section,
we first introduce a general form of the variational autoencoder (VAE) and
IWAE in the case of completely-observed data $\mathbf{X}$. Then, in Section
2.4.3, we apply the IWAE covariate structure to the deeply-learned GLM setting
and show how this representation naturally extends to the case where MCAR,
MAR, or MNAR missingness is observed in $\mathbf{X}$ when training deeply-
learned GLMs.

First, let $\mathbf{Z}$ be an $n\times d$ matrix, such that
$\mathbf{Z}=\\{\mathbf{z}_{1},\cdots,\mathbf{z}_{n}\\}$ and $\mathbf{z}_{i}$
is a latent vector of length $d$ pertaining to the $i^{th}$ sample latent
variable, and let $\mathbf{Z}$ represent a lower-dimensional representation or
subspace of $\mathbf{X}$. It is common practice to tune the value of $d$ as a
hyperparameter by choosing the optimal integer value that best fits the data,
as measured by some objective function. In a VAE, we assume
$\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ are i.i.d. samples from a multivariate
p.d.f or “generative model” $p_{\psi}(\mathbf{X}|\mathbf{Z})$ with
accompanying parameters $\psi$ that describes how $\mathbf{X}$ is generated
from the lower dimensional space $\mathbf{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 the training data. These
aspects are also aided through the use of embedded deep learning neural
networks, for example within $p_{\psi}(\mathbf{X}|\mathbf{Z})$, which also
facilitates its applicability to larger dimensions and complex datasets.

In a VAE with completely observed training data, one aims to maximize the
marginal log-likelihood as $\log p_{\psi}(\mathbf{X})=\log\int
p_{\psi}(\mathbf{X},\mathbf{Z})d\mathbf{Z}=\log\int
p_{\psi}(\mathbf{X}|\mathbf{Z})p(\mathbf{Z})d\mathbf{Z}$. However, this
quantity is also often intractable and difficult to maximize directly.
Therefore, VAE’s alternatively optimize an objective function called the
“Evidence Lower Bound” (ELBO), which lower bounds $\log p_{\psi}(\mathbf{X})$
and has the following form (Kingma & Welling, 2013):

$\displaystyle\mathcal{L}^{ELBO}(\theta,\psi)$
$\displaystyle=\operatorname{\mathbb{E}}_{\mathbf{Z}\sim
q_{\theta}(\mathbf{Z}|\mathbf{X})}\log\left[\frac{p_{\psi}(\mathbf{X}|\mathbf{Z})p(\mathbf{Z})}{q_{\theta}(\mathbf{Z}|\mathbf{X})}\right]$
(3) $\displaystyle\hat{\mathcal{L}}^{ELBO}_{K}(\theta,\psi)$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\log\left[\frac{p_{\psi}(\mathbf{X}|\tilde{\mathbf{Z}}_{k})p(\tilde{\mathbf{Z}}_{k})}{q_{\theta}(\tilde{\mathbf{Z}}_{k}|\mathbf{X})}\right].$
(4)

Here, $\mathcal{L}^{ELBO}(\theta,\psi)$ denotes the ELBO such that
$\mathcal{L}^{ELBO}(\theta,\psi)\leq\log p_{\psi}(\mathbf{X})$. Also let
$\hat{\mathcal{L}}^{ELBO}_{K}(\theta,\psi)$ denote the empirical approximation
to Eq. (3) computed by Monte Carlo integration, such that
$\mathcal{L}^{ELBO}(\theta,\psi)\approx\hat{\mathcal{L}}_{K}^{ELBO}(\theta,\psi)$
and $\tilde{\mathbf{Z}}_{1},\ldots,\tilde{\mathbf{Z}}_{K}$ are $K$ samples
drawn from $q_{\theta}(\mathbf{Z}|\mathbf{X})$, the variational approximation
of the true but intractable posterior $p_{\psi}(\mathbf{Z}|\mathbf{X})$, also
called the “recognition model”. Furthermore, denote $f_{\psi}(\mathbf{Z})$ and
$g_{\theta}(\mathbf{X})$ as the decoder and encoder feed forward neural
networks of the VAE, where $\psi$ and $\theta$ are the sets of weights and
biases pertaining to each of these neural networks, respectively. Given
$\mathbf{Z}$, $f_{\psi}(\mathbf{Z})$ outputs the distributional parameters
pertaining to $p_{\psi}(\mathbf{X}|\mathbf{Z})$.

In variational inference, $q_{\theta}(\mathbf{Z}|\mathbf{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_{\psi}(\mathbf{Z}|\mathbf{X})$. Variational inference is usually used in
tandem with amortization of the parameters where the neural network parameters
are shared across observations (Gershman & Goodman, 2014), allowing for
stochastic gradient descent (SGD) to be used for optimization of Eq. (4)
(Kingma & Welling, 2019). In practice, both
$q_{\theta}(\mathbf{Z}|\mathbf{X})$ and $p(\mathbf{Z})$ are typically assumed
to have simple forms, such as multivariate Gaussians with diagonal covariance
structures, and $q_{\theta}(\mathbf{Z}|\mathbf{X})$ is commonly assumed to be
factorizable, such that
$q_{\theta}(\mathbf{Z}|\mathbf{X})=\prod_{i=1}^{n}q_{\theta}(\mathbf{z}_{i}|\mathbf{x}_{i})$
(Kingma & Welling, 2019). Although one can specify a class of more complicated
distributions for $q_{\theta}(\mathbf{Z}|\mathbf{X})$ as long as they are
reparameterizable (Li et al., 2020; Strauss & Oliva, 2021, 2022), the
multivariate Gaussian with diagonal covariance structure is most often used,
following works by Burda et al. (2015) and Kingma & Welling (2013), due to the
convenience in sampling and computation (Kingma & Welling, 2019).

Let $(\hat{\theta}^{(t)},\hat{\psi}^{(t)})$ be the estimates of
$(\theta,\psi)$ 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 (Saxe et al., 2014; Murphy, 2016). Each
subsequent update $t\geq 1$ consists of two general steps to maximize
$\mathcal{L}(\theta,\psi)$. First, $K$ samples are drawn from
$q_{\hat{\theta}^{(t)}}(\mathbf{Z}|\mathbf{X})$ to compute the quantity in Eq.
(4), conditional on $\hat{\theta}^{(t)}$, similar to importance sampling.
Then, the so-called “reparametrization trick” is utilized to facilitate the
calculation of gradients of this approximation to obtain
$(\hat{\theta}^{t+1},\hat{\psi}^{t+1})$ using stochastic gradient descent
(Kingma & Welling, 2013). The networks $f_{\psi}(\mathbf{Z})$ and
$g_{\theta}(\mathbf{X})$ also allow the VAE to capture complex and non-linear
relationships between features in outputting the distributional parameters for
the generative and recognition models, respectively. This procedure may be
repeated for a fixed number of iterations, or may be terminated early due to
pre-specified convergence criteria (Prechelt, 1998). Kingma & Welling (2013)
provides additional details on the maximization procedure for VAEs.

The IWAE (Burda et al., 2015) is a generalization of the standard VAE. Both
the VAE and IWAE estimate $\log p_{\psi}(\mathbf{X})$ by drawing samples of
latent variables to estimate an expectation. However, while the VAE utilizes
$p_{\psi}(\mathbf{X},\mathbf{Z})/q_{\theta}(\mathbf{Z}|\mathbf{X})$ as the
importance weights in deriving the ELBO, the IWAE uses the average of $K$
importance weights in the integrand for a tighter lower bound of the marginal
log-likelihood (Burda et al., 2015). The resulting IWAE bound, corresponding
to the ELBO in Eq. (3), can be written as

$\displaystyle\mathcal{L}_{K}^{IWAE}(\theta,\psi)$
$\displaystyle=\operatorname{\mathbb{E}}_{\mathbf{Z}_{k}\sim
q_{\theta}(\mathbf{Z}|\mathbf{X})}\log\left[\frac{1}{K}\sum_{k=1}^{K}\frac{p_{\psi}(\mathbf{X}|\mathbf{Z}_{k})p(\mathbf{Z}_{k})}{q_{\theta}(\mathbf{Z}_{k}|\mathbf{X})}\right]$
(5) $\displaystyle\hat{\mathcal{L}}_{K}^{IWAE}(\theta,\psi)$
$\displaystyle=\log\left[\frac{1}{K}\sum_{k=1}^{K}\frac{p_{\psi}(\mathbf{X}|\tilde{\mathbf{Z}}_{k})p(\tilde{\mathbf{Z}}_{k})}{q_{\theta}(\tilde{\mathbf{Z}}_{k}|\mathbf{X})}\right].$
(6)

Importantly, although $K$ samples are drawn from $q(\mathbf{Z}|\mathbf{X})$ to
estimate the lower bound for both the VAE and IWAE, a VAE assumes a single
latent variable $\mathbf{Z}$ that is sampled $K$ times, wheras an IWAE assumes
$\mathbf{Z}_{1},\ldots,\mathbf{Z}_{K}$ are independent and identically
distributed (i.i.d.) latent variables, and each variable is sampled once from
$q(\mathbf{Z}|\mathbf{X})$. Typically, just one sample is drawn for each
latent variable to estimate the ELBO and IWAE bound. If $K=1$,
$\mathcal{L}_{1}^{IWAE}=\mathcal{L}^{VAE}$, and the IWAE corresponds exactly
to the standard VAE. For $K>1$, Burda et al. (2015) showed that $\log
p(\mathbf{X})\geq\hat{\mathcal{L}}_{K+1}^{IWAE}\geq\hat{\mathcal{L}}_{K}^{IWAE}$,
such that $\hat{\mathcal{L}}_{K}^{IWAE}\rightarrow\log p(\mathbf{X})$ as
$K\rightarrow\infty$ if
$p_{\psi}(\mathbf{X},\mathbf{Z})/q_{\theta}(\mathbf{Z}|\mathbf{X})$ is
bounded. Thus, the IWAE bound more closely approximates the true marginal log
likelihood when $K>1$ (Cremer et al., 2017), but the computational burden is
increased due to the increased number of samples. A visualization of the
workflow for an IWAE can be found in Appendix A3 of the supplementary
materials.

####  2.4.3 dlglm: Modeling X in the presence of missingness

Now, we extend the above framework to the deeply-learned GLM framework, where
features within $\mathbf{X}$ are partially observed during training. We
formally introduce the dlglm model to handle MNAR missingness in the context
of deeply-learned GLMs, as well as a variant of dlglm to specifically handle
MCAR and MAR missingness.

Let us define $q_{\theta}(\mathbf{Z},\mathbf{X}^{m})$ as the variational joint
posterior pertaining to $(\mathbf{Z},\mathbf{X}^{m})$. Then, we can factor
this variational joint posterior as
$q_{\theta}(\mathbf{Z},\mathbf{X}^{m})=q_{\theta_{1}}(\mathbf{Z}|\mathbf{X}^{o})q_{\theta_{2}}(\mathbf{X}^{m}|\mathbf{Z},\mathbf{X}^{o},\mathbf{R},\mathbf{Y})$.
Here, for $k=1,\ldots,K$, we assume $\mathbf{Z}_{k}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}q_{\theta_{1}}(\mathbf{Z}|\mathbf{X}^{o})$ similar to an IWAE,
and additionally assume $\mathbf{X}_{k}^{m}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}q_{\theta_{2}}(\mathbf{X}^{m}|\mathbf{Z},\mathbf{X}^{o},\mathbf{R},\mathbf{Y})$,
where each $\mathbf{X}_{k}^{m}$ has dimensionality $p_{miss}^{X}$. Here, we
assume that Y is generated by X, and thus it is redundant to utilize Y in the
part of the variational joint posterior pertaining to Z. Empirically, we
observed that including Y in the conditional, such that
$q(\mathbf{Z}|\mathbf{X}^{o},\mathbf{Y})$, did not have a significant impact.
Additionally, we note that the form
$q_{\theta_{2}}(\mathbf{X}^{m}|\mathbf{Z},\mathbf{X}^{o},\mathbf{R},\mathbf{Y})$
includes $\mathbf{Y}$, allowing for more accurate imputation of missing
values; however, we remove this term in the conditional in the context of
prediction, in order to predict $\mathbf{Y}$ in an unbiased manner.

We then utilize the class of factored variational posteriors, such that
$q_{\theta}(\mathbf{Z},\mathbf{X}^{m})=\prod_{i=1}^{n}q_{\theta}(\mathbf{z}_{i},\mathbf{x}_{i}^{m})$
and
$q_{\theta}(\mathbf{z}_{i},\mathbf{x}_{i}^{m})=q_{\theta_{1}}(\mathbf{z}_{i}|\mathbf{x}_{i}^{o})q_{\theta_{2}}(\mathbf{x}_{i}^{m}|\mathbf{z}_{i},\mathbf{x}_{i}^{o},\mathbf{r}_{i}^{X})$,
with $\theta=\\{\theta_{1},\theta_{2}\\}$. Then, denoting $\mathbf{z}_{ik}$
and $\mathbf{x}_{ik}^{m}$ as the $i^{th}$ observation vectors of
$\mathbf{Z}_{k}$ and $\mathbf{X}_{k}^{m}$, respectively, we have
$\mathbf{z}_{i1},\ldots,\mathbf{z}_{iK}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}q_{\theta_{1}}(\mathbf{z}_{i}|\mathbf{x}_{i}^{o})$ and
$\mathbf{x}_{i1}^{m},\ldots,\mathbf{x}_{iK}^{m}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}q_{\theta_{2}}(\mathbf{x}_{i}^{m}|\mathbf{z}_{i},\mathbf{x}_{i}^{o},\mathbf{r}_{i}^{X})$.
In this case, the lower bound, which we call the “dlglm bound”, can be derived
as follows:

$\displaystyle\log p$
${}_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{X}^{o},\mathbf{Y},\mathbf{R}^{X})=\sum_{i=1}^{n}\log
p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{x}_{i}^{o},y_{i},\mathbf{r}_{i}^{X})$
$\displaystyle=\sum_{i=1}^{n}\log\left[\iint
p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{x}_{i}^{o},\mathbf{x}_{i}^{m},y_{i},\mathbf{r}_{i}^{X},\mathbf{z}_{i})d\mathbf{z}_{i}d\mathbf{x}_{i}^{m}\right]$
$\displaystyle=\sum_{i=1}^{n}\log\operatorname{\mathbb{E}}_{(\mathbf{z}_{ik},\mathbf{x}_{ik}^{m})\sim
q_{\theta}(\mathbf{z}_{i},\mathbf{x}_{i}^{m})}\left[\frac{1}{K}\sum_{k=1}^{K}\frac{p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{x}_{i}^{o},\mathbf{x}_{ik}^{m},y_{i},\mathbf{r}_{i}^{X},\mathbf{z}_{ik})}{q_{\theta}(\mathbf{z}_{ik},\mathbf{x}_{ik}^{m})}\right]$
$\displaystyle\geq\sum_{i=1}^{n}\operatorname{\mathbb{E}}_{(\mathbf{z}_{ik},\mathbf{x}_{ik}^{m})\sim
q_{\theta}(\mathbf{z}_{i},\mathbf{x}_{i}^{m})}\log{\left[\frac{1}{K}\sum_{k=1}^{K}\frac{p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{x}_{i}^{o},\mathbf{x}_{ik}^{m},y_{i},\mathbf{r}_{i}^{X},\mathbf{z}_{ik})}{q_{\theta}(\mathbf{z}_{ik},\mathbf{x}_{ik}^{m})}\right]}=\mathcal{L}_{K}^{dlglm},$
(7)

Here, $\\{\psi,\boldsymbol{\beta},\pi,\phi,\theta\\}$ are the weights and
biases associated with the neural networks that output the parameters of the
distributions that are involved, $\alpha$ is the dispersion parameter
associated with the variance function of $\mathbf{Y}$, and
$\tilde{\mathbf{z}}_{ik}$ and $\tilde{\mathbf{x}}_{ik}^{m}$ are the samples
drawn from $q_{\theta_{1}}(\mathbf{z}_{i}|\mathbf{x}_{i}^{o})$, and
$q_{\theta_{2}}(\mathbf{x}_{i}^{m}|\mathbf{z}_{i},\mathbf{x}_{i}^{o},\mathbf{r}_{i}^{X})$,
respectively.

As discussed in Section 2.3, we use the selection model factorization of the
complete data log-likelihood, such that
$p_{\alpha,\boldsymbol{\beta},\pi,\psi,\phi}(\mathbf{x}_{i}^{o},\mathbf{x}_{i}^{m},y_{i},\mathbf{r}_{i}^{X},\mathbf{z}_{i})=p_{\alpha,\boldsymbol{\beta},\pi}(y_{i}|\mathbf{x}_{i})p_{\psi}(\mathbf{x}_{i}|\mathbf{z}_{i})p(\mathbf{z}_{i})p_{\phi}(\mathbf{r}_{i}^{X}|\mathbf{x}_{i},y_{i}).$
As before, we can remove $y_{i}$ from
$p_{\phi}(\mathbf{r}_{i}^{X}|\mathbf{x}_{i},y_{i})$ for unbiased prediction.
Then, applying this factorization to (7), we obtain the form of the estimate
of the “dlglm bound”, where the integral is estimated via Monte Carlo
integration:

$\hat{\mathcal{L}}_{K}^{dlglm}=\sum_{i=1}^{n}\log{\left[\frac{1}{K}\sum_{k=1}^{K}\frac{p_{\alpha,\boldsymbol{\beta},\pi}(y_{i}|\mathbf{x}_{i}^{o},\tilde{\mathbf{x}}_{ik}^{m})p_{\psi}(\mathbf{x}_{i}|\tilde{\mathbf{z}}_{ik})p(\tilde{\mathbf{z}}_{ik})p_{\phi}(\mathbf{r}_{i}^{X}|\mathbf{x}_{i}^{o},\tilde{\mathbf{x}}_{ik}^{m})}{q_{\theta_{1}}(\tilde{\mathbf{z}}_{ik}|\mathbf{x}_{i}^{o})q_{\theta_{2}}(\tilde{\mathbf{x}}_{ik}^{m}|\tilde{\mathbf{z}}_{ik},\mathbf{x}_{i}^{o},\mathbf{r}_{i}^{X})}\right]},$
(8)

Figure 2: Architecture of proposed dlglm method (Case x). Dark colored nodes
($X^{o},X^{m}=0,R^{X}$) represent deterministic values, lightly colored nodes
($Z^{\prime},X^{o\prime},X^{m\prime},R^{X\prime},Y^{\prime}$) represent
learned distributional parameters, and outlined (in red) nodes represent
sampled values. Orange cells correspond to latent variables $\mathbf{Z}$ and
$\mathbf{X}^{m}$. $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{K}$ and
$\mathbf{X}_{1}^{m},\ldots,\mathbf{X}_{K}^{m}$ are sampled from their
respective variational posteriors $q_{\theta_{1}}(\mathbf{Z}|\mathbf{X}^{o})$
and $q_{\theta_{2}}(\mathbf{X}^{m}|\mathbf{Z},\mathbf{R}^{X},\mathbf{X}^{o})$.

We see that this quantity closely resembles the lower bound of an IWAE, and,
similar to traditional VAEs, we utilize neural networks
$f_{\psi}(\mathbf{z}_{i})$, $g_{\theta_{1}}(\mathbf{x}_{i}^{o})$,
$g_{\theta_{2}}(\mathbf{z}_{i},\mathbf{x}_{i}^{o},\mathbf{r}_{i}^{X})$,
$s_{\boldsymbol{\beta},\pi}(\mathbf{x}_{i})$, and $h_{\phi}(\mathbf{x}_{i})$
to learn the values of the parameters of
$p_{\psi}(\mathbf{x}_{i}|\mathbf{z}_{i})$,
$q_{\theta_{1}}(\mathbf{z}_{i}|\mathbf{x}_{i}^{o})$,
$q_{\theta_{2}}(\mathbf{x}_{i}^{m}|\mathbf{z}_{i},\mathbf{x}_{i}^{o},\mathbf{r}_{i}^{X})$,
$p_{\alpha,\boldsymbol{\beta},\pi}(y_{i}|\mathbf{x}_{i})$, and
$p_{\phi}(\mathbf{r}_{i}^{X}|\mathbf{x}_{i})$. The associated weights and
biases of the neural networks $\\{\boldsymbol{\beta},\pi,\psi,\phi\\}$, as
well as the dispersion parameter $\alpha$ pertaining to
$p_{\alpha,\boldsymbol{\beta},\pi}(\mathbf{Y}|\mathbf{X})$ are updated using
stochastic gradient descent via the ADAM optimizer (Kingma & Ba, 2014).
Importantly, we call the neural network denoted by $h_{\phi}(\mathbf{x}_{i})$
the “missingness network”. The inclusion of this network allows us to learn a
model for the missingness mechanism, which is essential for accurate analysis
in the presence of MNAR or non-ignorable missingness. The architecture of
dlglm can be found in Figure 2. A pseudo-algorithm of dlglm can be found in
Appendix A4 of the supplementary materials. We limited our discussion in this
paper to Case x, where missingness exists only in $\mathbf{X}$ but not in Y;
however, the lower bound for dlglm can similarly be derived for the more
general Case xy as well, and this derivation can be found in Appendix A5 of
the supplementary materials.

We can obtain a variant of this method, which we call ignorably-missing dlglm
(idlglm), by assuming independence between $\mathbf{X}^{m}$ and $\mathbf{R}$
by omitting $\mathbf{r}_{i}^{X}$ from Equation 7, and removing
$p_{\phi}(\mathbf{r}_{i}^{X}|\mathbf{x}_{i}^{o},\tilde{\mathbf{x}}_{ik}^{m})$
and letting
$p_{\phi}(\tilde{\mathbf{x}}_{ik}^{m}|\tilde{\mathbf{z}}_{ik},\mathbf{x}_{i}^{o},\mathbf{r}_{i}^{X})\rightarrow
p_{\phi}(\tilde{\mathbf{x}}_{ik}^{m}|\tilde{\mathbf{z}}_{ik},\mathbf{x}_{i}^{o})$
in Equation 8. Whereas dlglm is better suited to handle MNAR, idlglm may be
more appropriate for the MCAR or MAR settings, where a missingness model need
not be specified.

In this paper, we are primarily interested in supervised learning. However,
following training, dlglm and idlglm can also perform imputation as in the
unsupervised learning architecture for handling missingness proposed by Lim et
al. (2021), although such imputation is not necessary for training,
coefficient estimation, or prediction. The single imputation procedure, and
additional computational details of dlglm and idlglm can be found in Appendix
A6 and A1 of the supplementary materials, respectively.

A recently published method by Ipsen et al. (2021) performs unsupervised
learning by similarly learning a missingness model in their neural network
framework to handle MNAR missingness. However, they assume that
$p(\mathbf{X}^{m}|\mathbf{Z})\approx
q(\mathbf{X}^{m}|\mathbf{Z},\mathbf{X}^{o},\mathbf{R})$. This may be an
oversimplification, as in the MNAR case, R cannot be assumed to be independent
of $\mathbf{X}^{m}$. Recent work by Ma & Zhang (2021) similarly performs
unsupervised learning under MNAR missingness, including an auxiliary fully-
observed variable to guarantee identifiability. However, they also make the
same simplifying assumption as Ipsen et al. (2021), which may not hold in
MNAR. In addition, both methods are designed for imputation, rather than
supervised tasks, and extending these methods for such tasks may be
nontrivial, especially for computationally intensive models. Issues of
identifiability in missing data applications often lead to issues of
convergence during model training (Beesley et al., 2019). We note that
although deriving the identifiability of dlglm is not focal point of this
paper, we consistently observed convergence in training the dlglm architecture
in various simulations and real data settings.

##  3 Numerical Examples

In this section, we evaluate the performance of dlglm and idlglm to analyze
each method’s performance in imputation, coefficient estimation, and
prediction tasks on simulated datasets under MCAR, MAR, and MNAR missingness
in Section 3.1. We also compare our methods to two commonly used approaches
for modeling missing data in the supervised setting, mean imputation and the
mice method for multiple imputation (Van Buuren & Groothuis-Oudshoorn, 2011).
We also compared performance in simulated data with two deep learning methods
that were recently published miwae (Mattei & Frellsen, 2019) and notmiwae
(Ipsen et al., 2021). To account for potential non-linearity and complex
relationships between features, in Section 3.2, we mask completely-observed
datasets obtained from the UCI Machine Learning Repository with varying
mechanisms of missingness on the predictors. Finally, in Section 3.3, we
perform prediction on the Bank Marketing dataset, which inherently contains
missingness in the predictors.

In all simulated and real data analyses, we tuned a variety of hyperparameters
for deep learning methods, including the number of hidden layers, the
dimensionality of the latent variable $\mathbf{Z}$, and the number of nodes
per hidden layer. For dlglm, we additionally tuned the number of hidden layers
in the missingness network separately, allowing the network to accurately
capture potentially complex nonlinear relationships in the missingness model.

A grid-search approach was used for training based upon discrete pre-specified
values, selecting the optimal combination of hyperparameters using the lower
bound computed on a held out validation set. The selected hyperparameters for
the simulated datasets, as well as the UCI and Bank Marketing datasets are
listed in Appendix B1 of the supplementary materials.

###  3.1 Simulated Data

####  3.1.1 Simulation Setup

We first utilized completely synthetic data to evaluate the performance of
each. Here, $\mathbf{X}$ is generated such that
$\mathbf{X}=normalize(\mathbf{Z}\mathbf{W}+\mathbf{B})+B_{0}$, where
$normalize(\cdot)$ takes an input matrix and standardizes each column to mean
0 and standard deviation 1, and $\mathbf{W}$ and $\mathbf{B}$ and are matrices
of dimensions $d\times p$ and $n\times p$, respectively, $\mathbf{Z}\sim
N_{d}(\mathbf{0},\mathbf{I})$, and $W_{lj}\sim N(0,0.5)$ and $B_{ij}\sim
N(0,1)$ for $i=1,\ldots,n$, $p=1,\ldots,p$, and $l=1,\ldots,d$, and $B_{0}=2$
is fixed. We also generated a binary response variable $\mathbf{Y}$ such that
logit$[$Pr$(\mathbf{Y}=1|\mathbf{X})]=\beta_{0}+\boldsymbol{\beta}\mathbf{X}$,
where $\boldsymbol{\beta}$ are drawn randomly from
$\\{-\frac{1}{4},\frac{1}{4}\\}$, and $\beta_{0}$ is chosen such that
approximately half of the sample are in either class. Values of $\mathbf{Y}$
are drawn from Bernoulli(Pr$(\mathbf{Y}=1|\mathbf{X})$).

We then simulate the missingness mask matrix $\mathbf{R}^{X}$ such that 50% of
features in $\mathbf{X}$ are partially observed, and 30% of the observations
for each of these features are missing. We generate $r_{ij}$ from the
Bernoulli distribution with probability equal to
$p(r_{ij_{m}}=1|\mathbf{x}_{i},y_{i},\boldsymbol{\phi})$, such that
$\text{logit}[p(r_{ij_{m}}=1|\mathbf{x}_{i},y_{i},\boldsymbol{\phi})]=\phi_{0}+\phi_{1}y_{i}+\boldsymbol{\phi}_{2}\mathbf{x}_{i}^{o}+\boldsymbol{\phi}_{3}\mathbf{x}_{i}^{m}$,
where $j_{m}=1,\ldots,p_{miss}^{X}$ index the missing features, $\phi_{1}$ is
the coefficient pertaining to the response,
$\boldsymbol{\phi}_{2}=\\{\phi_{21},\ldots,\phi_{2,p_{obs}^{X}}\\}$ are the
coefficients pertaining to the observed features, and
$\boldsymbol{\phi}_{3}=\\{\phi_{31},\ldots,\phi_{3,p_{miss}^{X}}\\}$ are those
pertaining to the missing features, where $p_{obs}^{X}$ and $p_{miss}^{X}$ are
the total number of features that are observed and missing, respectively, with
$p_{miss}^{X}=floor(0.5*p)$ and $p_{obs}^{X}=p-p_{miss}^{X}$. Here, we fixed
$\phi_{1}=0$, and drew nonzero values of
$\\{\boldsymbol{\phi}_{2},\boldsymbol{\phi}_{3}\\}$ from the log-normal
distribution with mean $\mu_{\phi}=5$, with $\log$ standard deviation
$\sigma_{\phi}=0.2$.

To evaluate the impact of the misspecification of the missingness mechanism on
model performance, $r_{ij_{m}}$ was simulated under each mechanism as follows:
(1) MCAR: $\\{\phi_{1},\boldsymbol{\phi}_{2},\boldsymbol{\phi}_{3}\\}=0$ (2)
MAR: Same as MCAR except $\phi_{2j_{o}}\neq 0$ for one completely-observed
feature $j_{o}$ (3) MNAR: Same as MCAR except $\phi_{3j_{m}}\neq 0$ for one
missing feature $j_{m}$. In this way, for each MAR or MNAR feature, the
missingness is dependent on just one feature. In each case, we used $\phi_{0}$
to control for an expected rate of missingness of $30\%$ in each partially-
observed feature. We note that for each these simulations, we utilize all
features in $\mathbf{X}$ as well as the response $\mathbf{Y}$ as input into
dlglm’s missingness network, although only one feature is involved under the
true missingness model. Additionally, we searched for the optimal variational
distributions of $q_{\theta_{1}}(\mathbf{Z}|\mathbf{X}^{o})$ and
$q_{\theta_{2}}(\mathbf{X}^{m}|\mathbf{Z},\mathbf{X}^{o},\mathbf{R})$ from a
class of Gaussian distributions with diagonal covariance structures, as
discussed in Section 2.4.2. We fixed $K=5$ during training, and increased $K$
to $500$ at test time.

We vary $n$ and $d$ such that $n=\\{10,000,100,000\\}$ and $d=\\{2,8\\}$, and
fix $p=50$. We simulated 5 datasets per simulation condition, spanning various
missingness mechanisms and values for $\\{n,d\\}$. We fix the values of
$\boldsymbol{\beta}$ at $0.25$ for each feature, and adjusted $\beta_{0}$ to
ensure equal proportions for the binary class response $\mathbf{Y}$. For each
simulation case, we partitioned the data into training, validation, and test
sets with ratio 8:1:1. For mice imputation, we averaged across 500 multiply-
imputed datasets to obtain a single imputed dataset. We note that we generated
$\mathbf{Y}$ by a linear transformation of $\mathbf{X}$ in these simulations
in order to facilitate fair comparisons with mice, which cannot account for
non-linear relationships between the features and the response. Because no
hyperparameter tuning is required, the validation set is not utilized for mice
and mean imputation.

We measured the performance of each method with respect to three different
tasks: imputation of missing values, coefficient estimation, and prediction.
Imputation performance was measured with respect to the truth on a single
imputed dataset by mean, dlglm and idlglm imputation, and on an average of
multiply-imputed datasets by mice. Coefficient estimation for mean, miwae,
notmiwae, and mice were based on downstream fitted GLM(s) on these imputed
dataset(s), where estimates were pooled using Rubin’s rules (Rubin, 2004) for
mice. For dlglm and idlglm, we estimated the coefficients by the weights and
bias $\boldsymbol{\beta}$ of the last layer of the
$s_{\boldsymbol{\beta},\pi}(\cdot)$ trained neural network. Here, we fixed the
number of hidden layers in $s_{\boldsymbol{\beta},\pi}(\cdot)$ to 0 to allow
for direct comparison with the other methods. A more complex prediction model
via a neural network can be learned by simply incorporating additional hidden
layers in $s_{\boldsymbol{\beta},\pi}(\cdot)$. We note that dlglm and idlglm
can estimate $\boldsymbol{\beta}$ without having to perform multiple
imputation and downstream modelling unlike mice, where fitting complex methods
such as neural networks each of the multiply-imputed datasets separately may
be computationally prohibitive.

After obtaining the coefficient estimates and trained models, we performed
prediction on the test set in two ways: 1) using the incomplete (predI) test
set, where the true values of $\mathbf{X}^{m}$ are not known at prediction
time, and 2) using the complete (predC) test set, where the true simulated
values of $\mathbf{X}^{m}$ are known at prediction time. These two ways
reflect the two realistic cases in which (1) missingness is present during
training time but complete data is available at prediction time, and (2)
missingness is present during both training and prediction time. For predI,
miwae, notmiwae, mice and mean imputation require an additional imputation
step on the test set before predicting $\mathbf{Y}$; for dlglm and idlglm, we
simply input the incomplete test set into the trained model without needing to
separately impute the test set, and we predict using the trained model. That
is, miwae, notmiwae, mice and mean imputation cannot generalize the trained
model to impute the test set, dlglm and idlglm provide a seamless framework to
utilize the already-trained model to impute and predict on a held-out test
set. For predC, we use the underlying true values of $\mathbf{X}^{m}$ to
predict on the test dataset.

Imputation error was measured by the average L1 distance between true and
imputed masked values in $\mathbf{X}$. Letting $\hat{\mathbf{X}}^{m}$ denote
the imputed masked values of the true $\mathbf{X}^{m}$ values of the missing
entries, we denote the average L1 distance is simply
$\frac{\mid\hat{\mathbf{X}}^{m}-\mathbf{X}^{m}\mid}{N_{miss}},$ where
$N_{miss}$ is the total number of missing entries in the dataset. Performance
in coefficient estimation was measured by the average percent bias (PB) of the
coefficient estimates compared to the truth, averaged across the $p$ features,
i.e.
$PB=100\times\frac{1}{p}\sum_{j=1}^{p}\frac{|\beta_{j}-\hat{\beta}_{j}|}{|\beta_{j}|}.$
Finally, predC and predI prediction error was measured by the average L1
distance between predicted and true values of the probabilities of class
membership Pr$(\mathbf{Y}=1|\mathbf{X})$ in the test set.

In order to assess the sensitivity of the performance of these methods to the
specification of the missingness model used to synthetically mask the data, we
also repeated the analyses on data with missingness mask simulated by the
following:
$\text{logit}[p(r_{ij_{m}}=1|\mathbf{x}_{i},y_{i},\boldsymbol{\phi})]=\phi_{0}+\boldsymbol{\phi}_{2}\text{log}(\mathbf{x}_{i}^{o}+min(\mathbf{x}_{i}^{o}))+\boldsymbol{\phi}_{3}\text{log}(\mathbf{x}_{i}^{m}+min(\mathbf{x}_{i}^{m})),$
such that for the MAR and MNAR missingness cases, the missingness was
dependent on the log of one of the completely or partially observed features.
We denote this set of simulation conditions the “nonlinear missingness” case,
where the missingness was simulated from the specified nonlinear logistic
regression model. We show the results of this analysis in Appendix B2 of the
supplementary materials.

####  3.1.2 Simulation Results

Figure 3: Simulation results with $n=10,000$ and $p=50$, varying $d=2$ (top 4)
and $d=8$ (bottom 4). In each quadrant, we measure imputation accuracy by the
average L1 distance between imputed vs true values in $\mathbf{X}$ (top-left),
coefficient estimation accuracy by the average percent bias (PB) of the
estimates $\hat{\boldsymbol{\beta}}$ compared to the truth (top-right), and
prediction accuracy by the average L1 distance between the predicted and true
probabilities of class 1 membership of $\mathbf{Y}$ using the true unmasked
test set (predC, bottom-left) and the incomplete test set (predI, bottom-
right). In predI, we first impute missing values of the test data for mean,
miwae, notmiwae, and mice imputation, and we input the incomplete test set as-
is for dlglm and idlglm. Figure 4: Simulation results with $n=100,000$ and
$p=50$, varying $d=2$ (top 4) and $d=8$ (bottom 4). In each quadrant, we
measure imputation accuracy by the average L1 distance between imputed vs true
values in $\mathbf{X}$ (top-left), coefficient estimation accuracy by the
average percent bias (PB) of the estimates $\hat{\boldsymbol{\beta}}$ compared
to the truth (top-right), and prediction accuracy by the average L1 distance
between the predicted and true probabilities of class 1 membership of
$\mathbf{Y}$ using the true unmasked test set (predC, bottom-left) and the
incomplete test set (predI, bottom-right). In predI, we first impute missing
values of the test data for mean, miwae, notmiwae, and mice imputation, and we
input the incomplete test set as is for dlglm and idlglm.

Figures 3 and 4 illustrate the simulation results pertaining to imputation
accuracy, coefficient estimation, and prediction accuracy for the condition
$p=50$. We see that across all combinations of $\\{n,d\\}$ and mechanisms of
missingness, mean imputation consistently performs poorly in imputation,
coefficient estimation, and prediction, while mice and idlglm perform
comparably in these metrics. Also, we note that under MNAR missingness, dlglm
generally yields the lowest imputation and prediction error, as well as
percent bias across all simulation cases. Under MAR missingness, dlglm
performs comparably to idlglm and mice. This shows the ability of dlglm to
learn an accurate model of the missingness, even under severe
overparametrization of the missingness model (model need not be specified for
ignorable missingness). However, due to the complexity of the model, we see
that dlglm does generally perform poorly compared to idlglm and mice under
MCAR missingness, when $n=10,000$, although it still performs comparably to
other methods when the sample size is very large ($n=100,000$). As one may
expect, prediction performance using the incomplete data (predI) was poorer
than prediction performance using the complete data (predC) for all methods.

We additionally show results pertaining to $p=25$ in Appendix B2 of the
supplementary materials. We similarly found that dlglm performed best under
MNAR missingness, and comparably to idlglm, mice, miwae, and notmiwae under
MCAR and MAR missingness.

###  3.2 Real Data with Simulated Missingness

Next, we analyzed 3 completely-observed, large datasets from the UCI Machine
Learning Repository (Dua & Graff, 2017) that contained a specific response
variable of interest, in order to preserve non-linearity and interactions
between observed features. Unlike the simulated datasets, these UCI datasets
don’t follow a specific distribution that may be leveraged to inform a
supervised learning method. The DRYBEAN dataset contains 16 features
describing 13,611 images of dry beans taken with a high-resolution camera, and
the response variable of interest was the type of dry bean each image
represents, with 7 different possible types of beans. The LETTER dataset
contains 16 attributes of 20,000 black-and-white pixel images, each displaying
one English letter (A to Z). Finally, the SHUTTLE dataset contains 9 numerical
attributes pertaining to 58,000 shuttle stat logs (observations), which are
classified into 7 different categories. Due to a low sample size in 4 of the 7
categories, we pre-filtered the observations pertaining to these categories
out of the dataset, and the resulting dataset contained 57,756 observations of
3 categories. In each of these datasets, the response variable was categorical
with greater than two levels. Additional information regarding these datasets,
and how to obtain the raw data files can be found in Appendix C of the
supplementary materials.

We then simulated the missingness mask $\mathbf{R}^{X}$ with MCAR, MAR, and
MNAR patterns of missingness in the manner described in Section 3.1.1. We
split the samples in each dataset by a similar 8:1:1 ratio of
training/validation/test samples. In the test set samples, we then imputed the
missing values and predicted the response variables with each method in a
manner similar to Section 3.1.1. For dlglm and idlglm, We account for
potential nonlinear relationships between the covariates and response by
allowing the number of hidden layers in $s_{\boldsymbol{\beta},\pi}(\cdot)$ to
be nonzero in hyperparameter tuning. Then, we compared imputation and
prediction accuracy on each dataset, under each mechanism of missingness.
Since the underlying true probabilities of class membership were unavailable,
we measured prediction accuracy by the Cohen’s kappa metric on the complete
(kappaC) and incomplete (kappaI) test set, with predicted class determined by
the maximum predicted probability of membership. This metric measures how
accurately a categorical variable was predicted, with a value of -1 indicating
worst possible performance, and a value of 1 indicating perfect concordance
with the truth.

Figure 5: Imputation (top row) and prediction results from predC (middle row)
and predI (bottom row) from comparative methods run on 3 large datasets from
the UCI Machine Learning Repository: DRYBEAN, LETTER, and SHUTTLE (columns,
left to right). Imputation error was measured by the average L1 distance
between true and imputed entries, with lower values indicating better
performance, and prediction performance was measured by the Cohen’s kappa
metric for both predC (kappaC) and predI (kappaI), with higher values
indicating better performance.

Results from the imputation and prediction analyses on these datasets can be
found in Figure 5. We found that, as in the simulations, mean imputation
performed most poorly in both imputation and downstream prediction, while
dlglm tended to perform best in the MNAR cases, and performed comparably to
mice and idlglm under the MCAR and MAR cases. This further validates our
claims under a more realistic setting, where the true data generation
mechanism may be unknown. We also see that under both MCAR and MAR
missingness, mice performed worse than idlglm in prediction on the LETTER and
SHUTTLE datasets. The mice model has been known to break down under nonlinear
relationships between the features (Van Buuren, 2018), as may be the case in
real-world datasets like the ones being examined. Using neural networks to
model the data generation process allows idlglm to better model potential
nonlinear relationships between features, allowing for more accurate
prediction.

Interestingly, all of the algorithms performed similarly in prediction on the
DRYBEAN dataset. We found that this dataset contained extremely high levels of
correlation between the variables (see Web Appendix C of the supplementary
materials). When features containing missingness are highly correlated to
other fully-observed features, such missingness may not truly reflect the MNAR
scenario (Hapfelmeier et al., 2012). This is because there exist fully-
observed features that are highly correlated with the missing features, and
ignorably-missing data methods like idlglm may gather information about the
missing entries from these correlated, fully-observed features without having
to explicitly model the mechanism of missingness. Still, idlglm and dlglm
imputed missing entries much more accurately than mean and mice under MAR and
MNAR. Interestingly, we also see that idlglm performed similarly to dlglm
under MNAR in this dataset.

We additionally performed similar analyses on 5 other smaller UCI datasets,
and these results can be found in Appendix D of the supplementary materials.
We found that under small sample size settings, performance via dlglm was more
variable under MNAR. We suggest use of dlglm when the data contains at least
10,000 samples, as the model may be too complex to be accurately trained under
smaller sample sizes.

###  3.3 Bank Marketing Dataset

Finally, we performed prediction on the Bank Marketing dataset from the UCI
Machine Learning Repository. This dataset contained 41,188 observations of 20
different attributes that were obtained based on direct phone calls from a
Portuguese banking institution as part of a promotion campaign for a term
deposit subscription (Moro et al., 2014). The response variable of interest
was a fully-observed binary measure of whether the client subscribed a term
deposit. Of the 20 attributes, we removed 1 attribute as directed from the
manual due to perfect correlation with the response variable, and removed 3
other attributes that were deemed irrelevant to the prediction task: month of
contact, day of contact, and communication type (cell phone or telephone).

Missingness was present in 8 of the 16 attributes: type of job, marital
status, level of education, whether the client had a credit in default,
whether the client had a housing loan, whether the client had a personal loan,
number of days since the client was contacted in a previous campaign, and
outcome of the previous campaign. The remaining 8 attributes were fully-
observed: age of client, number of contacts during this campaign, number of
contacts before this campaign, employment variation rate, consumer price
index, consumer confidence index, euribor 3 month rate, and employee number.
The global rate of missingness was about 13.3%. The response variable of
interest was collected by additional follow-up calls to confirm whether the
client subscribed to the product. Additional information regarding the bank
marketing dataset, and how to obtain the raw data files can be found in
Appendix E of the supplementary materials.

This type of dataset reflects the most realistic situation in practice, where
missingness exists in a dataset and one has no prior knowledge of either the
relationships between the features and the response, or the underlying
mechanism of the missingness. We divided the dataset into the 8:1:1 training,
validation, and test set ratio, and performed prediction as before. Because
neither the data nor the missingness was simulated, we compared just the predI
prediction performance across the methods.

In order to more deeply dive into this real data example, we assessed the
prediction performance for dlglm and idlglm in the context of prediction
(excluding $\mathbf{Y}$ from neural networks $g_{\theta_{2}}(\cdot)$ and
$h_{\phi}(\cdot)$) and imputation (including $\mathbf{Y}$, denoted by dlglmy
and idlglmy).

| AUC | PPV | kappa | F1
---|---|---|---|---
dlglm | 0.778 | 0.475 | 0.397 | 0.470
dlglmy | 0.880 | 0.481 | 0.445 | 0.516
idlglm | 0.791 | 0.475 | 0.407 | 0.460
idlglmy | 0.779 | 0.446 | 0.411 | 0.488
mean | 0.769 | 0.448 | 0.385 | 0.46
mice | 0.771 | 0.455 | 0.396 | 0.471
Table 1: Results from prediction analyses on the Bank Marketing dataset from
the UCI Machine Learning Repository. We measured concordance between the true
and predicted binary response by 4 metrics: Area Under ROC Curve (AUC),
Positive Predictivity (PPV), Cohen’s kappa (kappa), and F1 score. For dlglmy
and idlglmy, $\mathbf{Y}$ was included in the input for neural networks
$g_{\theta_{2}}(\cdot)$ and $h_{\phi}(\cdot)$, in order to allow for more
accurate imputed values.

Table 1 shows the results from these prediction analyses. We measured
prediction performance of the binary response variable by 4 metrics: Area
Under the ROC Curve (AUC), Positive Predictivity (PPV), Cohen’s kappa (kappa),
and the F1 metric. The formulas for PPV and F1 metrics are given in Appendix F
of the supplementary materials. For each metric, a larger value represents
greater concordance between the true and predicted response. We see that
although dlglmy yielded a significantly greater performance in prediction via
all metrics, dlglm does not significantly outperform idlglm. The similar
performance between dlglm and idlglm may indicate that the real mechanism of
missingness in this data may not be MNAR, although this claim is not testable
in practice.

Additionally, the trained dlglm model chose 0 hidden layers in the
$s_{\boldsymbol{\beta},\boldsymbol{\pi}}(\mathbf{X})$ neural network in the
optimal model, such that
$s_{\boldsymbol{\beta},\boldsymbol{\pi}}(\mathbf{X})\rightarrow
s_{\boldsymbol{\beta}}(\mathbf{X})$. Therefore, the weights of that neural
network $\boldsymbol{\beta}$ coincide exactly with the coefficient estimates
of a classic generalized linear model, i.e.
$\text{logit}[Pr(Y=1)]=\mathbf{X}\boldsymbol{\beta}$. The features in the
dataset with the largest effects on the probability of a client subscribing to
a term deposit were employment variation rate (0.538), age of client (0.508),
and whether the client had a personal loan (-0.402). Specifically, a client
was more likely to subscribe if the company experienced higher levels of
variation in employment and if the client were older, while a client was less
likely to subscribe if they had a personal loan.

##  4 Discussion

In this paper, we introduced a novel deep learning method called Deeply-
learned Generalized Linear Model with Missing Data (dlglm), which is able to
perform coefficient estimation and prediction in the presence of missing not
at random (MNAR) data. dlglm utilizes a deep learning neural network
architecture to model the generation of the data matrix $\mathbf{X}$, as well
as the relationships between the response variable $\mathbf{Y}$ and
$\mathbf{X}$ and between the missingness mask $\mathbf{R}$ and $\mathbf{X}$.
In this way, we are able to (1) generalize the traditional GLM to account for
complex nonlinear interactions between the features, and (2) account for
ignorable and non-ignorable forms of missingness in the data. We also showed
through simulations and real data analyses that dlglm performs better in
coefficient estimation and prediction in the presence of MNAR missingness than
other impute-then-regress methods, like mean and mice imputation. Furthermore,
we found that dlglm was generally robust to the mechanism of missingness,
performing comparably well to mice and idlglm under MCAR and MAR settings.
Still, it is recommended to utilize idlglm when assuming the missingness is
ignorable, given that the missingness model that is learned in dlglm is not
necessary in this setting.

Supervised learning algorithms such as dlglm and idlglm can be particularly
useful in analyzing real-life data in the presence of missingness. In reality,
the mechanism underlying missing values cannot be explicitly known or tested,
but dlglm may allow flexibility to evaluate multiple assumptions regarding the
missingness mechanism. Furthermore, whereas impute-then-regress methods may
typically require fully-observed observations at test time for prediction,
dlglm and idlglm can predict the response of interest using partially-observed
observations. This provides a convenient workflow, where a user need not
separately re-impute the prediction set at test time.

In this paper, we focused specifically on the case of univariate response
$\mathbf{Y}$. dlglm and idlglm can be generalized to the multivariate Y case
by (1) including $\mathbf{Y}$ in the existing IWAE structure and (2) expanding
the neural network $s_{\boldsymbol{\beta},\pi}(\mathbf{x}_{i})$ to account for
all $q$ responses in $\mathbf{Y}$, and utilizing samples of $\mathbf{Z}$ as
additional input into this network such that
$s_{\boldsymbol{\beta},\pi}(\mathbf{x}_{i})\rightarrow
s_{\boldsymbol{\beta},\pi}(\mathbf{x}_{i},\mathbf{z}_{i})$. By doing this, we
account for multivariate $\mathbf{Y}$, outputting additional parameters
pertaining to the newly-specified distribution of
$p_{\boldsymbol{\beta},\pi}(\mathbf{y}_{i}|\mathbf{x}_{i},\mathbf{z}_{i})$ and
modelling correlation of Y by the learned latent structure. We leave this as
an extension of our method.

SUPPLEMENTARY MATERIAL

Supplementary Materials:

Additional details of the dlglm algorithm and the datasets used in this paper.
(pdf)

R-package for dlglm:

R-package dlglm containing code to perform the diagnostic methods described in
the article. The package can be downloaded from
https://github.com/DavidKLim/dlglm (website)

R Paper repo for reproducibility:

Github repository to replicate all analyses from this paper can be found here:
https://github.com/DavidKLim/dlglm_Paper (website)

## References

  * (1)
  * Beesley et al. (2019) Beesley, L. J., Taylor, J. M. & Little, R. J. (2019), ‘Sequential imputation for models with latent variables assuming latent ignorability’, Australian & New Zealand Journal of Statistics 61(2), 213–233.
  * Bottou (2012) Bottou, L. (2012), Stochastic gradient descent tricks, in ‘Neural networks: Tricks of the trade’, Springer, pp. 421–436.
  * Burda et al. (2015) Burda, Y., Grosse, R. & Salakhutdinov, R. (2015), ‘Importance Weighted Autoencoders’, arXiv e-prints p. arXiv:1509.00519.
  * Chen et al. (2019) Chen, D., Liu, S., Kingsbury, P., Sohn, S., Storlie, C. B., Habermann, E. B., Naessens, J. M., Larson, D. W. & Liu, H. (2019), ‘Deep learning and alternative learning strategies for retrospective real-world clinical data’, npj Digital Medicine 2(1).
  * Cremer et al. (2017) Cremer, C., Morris, Q. & Duvenaud, D. (2017), ‘Reinterpreting Importance-Weighted Autoencoders’, arXiv e-prints p. arXiv:1704.02916.
  * Diggle & Kenward (1994) Diggle, P. & Kenward, M. G. (1994), ‘Informative drop-out in longitudinal data analysis’, Applied Statistics 43(1), 49.
  * Dormehl (2019) Dormehl, L. (2019), ‘What is an artificial neural network? here’s everything you need to know’.
https://www.digitaltrends.com/cool-tech/what-is-an-artificial-neural-network/

  * Dua & Graff (2017) Dua, D. & Graff, C. (2017), ‘UCI machine learning repository’.
http://archive.ics.uci.edu/ml

  * Friedman et al. (2010) Friedman, J., Hastie, T. & Tibshirani, R. (2010), ‘Regularization paths for generalized linear models via coordinate descent’, Journal of Statistical Software 33(1).
  * Gershman & Goodman (2014) Gershman, S. J. & Goodman, N. D. (2014), Amortized inference in probabilistic reasoning, in ‘CogSci’.
  * Ghorbani & Zou (2018) Ghorbani, A. & Zou, J. Y. (2018), Embedding for informative missingness: Deep learning with incomplete data, in ‘2018 56th Annual Allerton Conference on Communication, Control, and Computing (Allerton)’, IEEE, pp. 437–445.
  * Guo & Gelfand (1990) Guo, H. & Gelfand, S. B. (1990), Analysis of gradient descent learning algorithms for multilayer feedforward neural networks, in ‘29th IEEE Conference on Decision and Control’, IEEE, pp. 1751–1756.
  * Hapfelmeier et al. (2012) Hapfelmeier, A., Hothorn, T., Ulm, K. & Strobl, C. (2012), ‘A new variable importance measure for random forests with missing data’, Statistics and Computing 24(1), 21–34.
  * Holland & Welsch (1977) Holland, P. W. & Welsch, R. E. (1977), ‘Robust regression using iteratively reweighted least-squares’, Communications in Statistics - Theory and Methods 6(9), 813–827.
  * Hoogland et al. (2020) Hoogland, J., Barreveld, M., Debray, T. P. A., Reitsma, J. B., Verstraelen, T. E., Dijkgraaf, M. G. W. & Zwinderman, A. H. (2020), ‘Handling missing predictor values when validating and applying a prediction model to new patients’, Statistics in Medicine 39(25), 3591–3607.
  * Ibrahim et al. (2005) Ibrahim, J. G., Chen, M.-H., Lipsitz, S. R. & Herring, A. H. (2005), ‘Missing-data methods for generalized linear models’, Journal of the American Statistical Association 100(469), 332–346.
  * Ibrahim & Molenberghs (2009) Ibrahim, J. G. & Molenberghs, G. (2009), ‘Missing data methods in longitudinal studies: a review’, TEST 18(1), 1–43.
  * Ipsen et al. (2021) Ipsen, N. B., Mattei, P.-A. & Frellsen, J. (2021), How to deal with missing data in supervised deep learning?, in ‘International Conference on Learning Representations’.
  * Kingma & Ba (2014) Kingma, D. P. & Ba, J. (2014), ‘Adam: A Method for Stochastic Optimization’, arXiv e-prints p. arXiv:1412.6980.
  * Kingma & Welling (2013) Kingma, D. P. & Welling, M. (2013), ‘Auto-Encoding Variational Bayes’, arXiv e-prints p. arXiv:1312.6114.
  * Kingma & Welling (2019) Kingma, D. P. & Welling, M. (2019), ‘An Introduction to Variational Autoencoders’, arXiv e-prints p. arXiv:1906.02691.
  * Li et al. (2020) Li, Y., Akbar, S. & Oliva, J. (2020), Acflow: Flow models for arbitrary conditional likelihoods, in ‘International Conference on Machine Learning’, PMLR, pp. 5831–5841.
  * Lim et al. (2021) Lim, D. K., Rashid, N. U., Oliva, J. B. & Ibrahim, J. G. (2021), ‘Handling Non-ignorably Missing Features in Electronic Health Records Data Using Importance-Weighted Autoencoders’, arXiv e-prints p. arXiv:2101.07357.
  * Lipsitz & Ibrahim (1996) Lipsitz, S. R. & Ibrahim, J. G. (1996), ‘A conditional model for incomplete covariates in parametric regression models’, Biometrika 83(4), 916–922.
  * Little & Rubin (2002) Little, R. J. A. & Rubin, D. B. (2002), Statistical Analysis with Missing Data, John Wiley & Sons, Inc.
  * Lopez et al. (2018) Lopez, R., Regier, J., Cole, M. B., Jordan, M. I. & Yosef, N. (2018), ‘Deep generative modeling for single-cell transcriptomics’, Nature Methods 15(12), 1053–1058.
  * Ma & Zhang (2021) Ma, C. & Zhang, C. (2021), ‘Identifiable generative models for missing not at random data imputation’, Advances in Neural Information Processing Systems 34, 27645–27658.
  * Mattei & Frellsen (2019) Mattei, P.-A. & Frellsen, J. (2019), MIWAE: Deep generative modelling and imputation of incomplete data sets, in K. Chaudhuri & R. Salakhutdinov, eds, ‘Proceedings of the 36th International Conference on Machine Learning’, Vol. 97 of Proceedings of Machine Learning Research, PMLR, Long Beach, California, USA, pp. 4413–4423.
  * McCullagh & Nelder (2019) McCullagh, P. & Nelder, J. A. (2019), Generalized linear models, Routledge.
  * Moro et al. (2014) Moro, S., Cortez, P. & Rita, P. (2014), ‘A data-driven approach to predict the success of bank telemarketing’, Decision Support Systems 62, 22–31.
  * Murphy (2016) Murphy, J. (2016), ‘An overview of convolutional neural network architectures for deep learning’, Microway Inc pp. 1–22.
  * Nelder & Wedderburn (1972) Nelder, J. A. & Wedderburn, R. W. M. (1972), ‘Generalized linear models’, Journal of the Royal Statistical Society. Series A (General) 135(3), 370–384.
  * Prechelt (1998) Prechelt, L. (1998), Early stopping-but when?, in ‘Neural Networks: Tricks of the trade’, Springer, pp. 55–69.
  * Qi & Wu (2003) Qi, M. & Wu, Y. (2003), ‘Nonlinear prediction of exchange rates with monetary fundamentals’, Journal of Empirical Finance 10(5), 623–640.
  * Razzak et al. (2017) Razzak, M. I., Naz, S. & Zaib, A. (2017), Deep learning for medical image processing: Overview, challenges and the future, in ‘Lecture Notes in Computational Vision and Biomechanics’, Springer International Publishing, pp. 323–350.
  * Rubin (1976) Rubin, D. B. (1976), ‘Inference and missing data’, Biometrika 63(3), 581–592.
  * Rubin (2004) Rubin, D. B. (2004), Multiple imputation for nonresponse in surveys, Vol. 81, John Wiley & Sons.
  * Saxe et al. (2014) Saxe, A. M., Mcclelland, J. L. & Ganguli, S. (2014), Exact solutions to the nonlinear dynamics of learning in deep linear neural network, in ‘In International Conference on Learning Representations’.
  * Strauss & Oliva (2021) Strauss, R. & Oliva, J. B. (2021), ‘Arbitrary conditional distributions with energy’, Advances in Neural Information Processing Systems 34, 752–763.
  * Strauss & Oliva (2022) Strauss, R. & Oliva, J. B. (2022), ‘Posterior matching for arbitrary conditioning’, Advances in Neural Information Processing Systems 35, 18088–18099.
  * Stubbendick & Ibrahim (2003) Stubbendick, A. L. & Ibrahim, J. G. (2003), ‘Maximum likelihood methods for nonignorable missing responses and covariates in random effects models’, Biometrics 59(4), 1140–1150.
  * Svozil et al. (1997) Svozil, D., Kvasnicka, V. & Pospichal, J. (1997), ‘Introduction to multi-layer feed-forward neural networks’, Chemometrics and intelligent laboratory systems 39(1), 43–62.
  * Tran et al. (2019) Tran, M.-N., Nguyen, N., Nott, D. & Kohn, R. (2019), ‘Bayesian deep net GLM and GLMM’, Journal of Computational and Graphical Statistics 29(1), 97–113.
  * Van Buuren (2018) Van Buuren, S. (2018), Flexible imputation of missing data, CRC press.
  * Van Buuren & Groothuis-Oudshoorn (2011) Van Buuren, S. & Groothuis-Oudshoorn, K. (2011), ‘mice: Multivariate imputation by chained equations in r’, Journal of statistical software 45(1), 1–67.
  * Wells et al. (2013) Wells, B. J., Nowacki, A. S., Chagin, K. & Kattan, M. W. (2013), ‘Strategies for handling missing data in electronic health record derived data’, eGEMs (Generating Evidence &: Methods to improve patient outcomes) 1(3), 7.