# Semiparametric fractional imputation using Gaussian mixture models for
handling multivariate missing data

Hejian Sang
Google Inc
and
Jae Kwang Kim
Department of Statistics, Iowa State University,Ames  Acknowledgment This work
is done when the author was in Iowa State University

###### Abstract

Item nonresponse is frequently encountered in practice. Ignoring missing data
can lose efficiency and lead to misleading inference. Fractional imputation is
a frequentist approach of imputation for handling missing data. However, the
parametric fractional imputation of Kim (2011) may be subject to bias under
model misspecification. In this paper, we propose a novel semiparametric
fractional imputation method using Gaussian mixture models. The proposed
method is computationally efficient and leads to robust estimation. The
proposed method is further extended to incorporate the categorical auxiliary
information. The asymptotic model consistency and $\sqrt{n}$\- consistency of
the semiparametric fractional imputation estimator are also established. Some
simulation studies are presented to check the finite sample performance of the
proposed method.

Keywords: Item nonresponse; Robust estimation; Variance estimation.

##  1 Introduction

Missing data is frequently encountered in survey sampling, clinical trials and
many other areas. Imputation can be used to handle item nonresponse and
several imputation methods have been developed in the literature. Motivated
from a Bayesian perspective, Rubin (1996) proposed multiple imputation to
create multiple complete data sets. Alternatively, under the frequentist
framework, fractional imputation (Kim and Fuller, 2004; Kim, 2011) makes one
complete data with multiple imputed values and their corresponding fractional
weights. Little and Rubin (2002) and Kim and Shao (2013) provided
comprehensive overviews of the methods for handling missing data.

For multivariate missing data with arbitrary missing patterns, valid
imputation methods should preserve the correlation structure in the imputed
data. Judkins et al. (2007) proposed an iterative hot deck imputation
procedure that is closely related to the data augmentation algorithm of Tanner
and Wong (1987) but they did not provide variance estimation. Other imputation
procedures for multivariate missing data include the multiple imputation
approaches of Raghunathan et al. (2001) and Murray and Reiter (2016), and
parametric fractional imputation of Kim (2011). The approaches of Judkins et
al. (2007) and Raghunathan et al. (2001) are based on conditionally specified
models and the imputation from the conditionally specified model is
generically subject to the model compatibility problem (Chen, 2010; Liu et
al., 2013; Bartlett et al., 2015). Conditional models for the different
missing patterns calculated directly from the observed patterns may not be
compatible with each other. The parametric fractional imputation uses the
joint distribution to create imputed values and does not suffer from model
compatibility problems.

Note that parametric imputation requires correct model specification.
Nonparametric imputation methods, such as kernel regression imputation (Cheng,
1994; Wang and Chen, 2009), are robust but may be subject to the curse of
dimensionality. Hence, it is important, often critical, to develop a unified,
robust and efficient imputation method that can be used for general purpose
estimation. The proposed semiparametric method fills in this important gap by
considering a more flexible model for imputation.

In this paper, to achieve robustness against model misspecification, we
develop an imputation procedure based on Gaussian mixture models. Gaussian
mixture model is a very flexible model that can be used to handle outliers,
heterogeneity and skewness. McLachlan and Peel (2004) and Bacharoglou (2010)
argued that any continuous distribution can be approximated by a finite
Gaussian mixture distribution. The proposed method using Gaussian mixture
model makes a nice compromise between efficiency and robustness. It is
semiparametric in the sense that the number of mixture components is chosen
automatically from the data. The computation for parameter estimation in our
proposed method is based on EM algorithm and its implementation is relatively
simple and efficient.

We note that Di Zio et al. (2007) also proposed to use Gaussian mixture model
to impute missing data. However, variance estimation and choice of the mixture
component are not discussed in Di Zio et al. (2007). Elliott and Stettler
(2007) and Kim et al. (2014) introduced the multiple imputation using mixture
models. Instead of multiple imputation, we use fractional imputation for
general-purpose estimation. We provide a completely theoretical justification
for consistency of the proposed imputation method. The variance estimator and
the model selection for the number of mixture component are also carefully
discussed and demonstrated in numerical studies. The proposed method is
further extended to handle mixed type data including categorical variable in
Section 5. By allowing the proportion vector of mixture component to depend on
categorical auxiliary variable, the proposed fractional imputation using
Gaussian mixture models can incorporate the observed categorical variables and
provide a very flexible tool for imputation.

The paper is structured as follows. The setup of the problem is introduced and
a short review of fractional imputation are presented in Secction 2. In
Section 3, the proposed semiparametric method and its algorithm for
implementation are introduced. Some asymptotic results are presented in
Section 4. In Section 5, the proposed method is further extended to handle
mixed type data. Some numerical studies and a real data application are
presented to show the performance of the proposed method in Section 6 and
Section 7, respectively. In Section 8, some concluding remarks are made. The
technical derivations and proof are presented in Appendix.

##  2 Basic Setup

Consider a $p$-dimensional vector of study variable
$Y=(y_{1},y_{2},\cdots,y_{p})$. Suppose that
$\left\\{Y_{1},Y_{2},\cdots,Y_{n}\right\\}$ are $n$ independent and
identically distributed realizations of the random vector $Y$. In this paper,
we use the upper case to represent vector or matrix and the lower case to
denote the elements within vector or matrix. Assume that we are interested in
estimating parameter $\theta\in\Theta$, which is defined through
$E\left\\{U(\theta;Y)\right\\}=0$, where $U(\cdot;Y)$ is the estimating
function of $\theta$. With no missingness, a consistent estimator of $\theta$
can be obtained by the solution to

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}U(\theta;Y_{i})=0.$ (1)

To avoid unnecessary details, we assume that the solution to (1) exists
uniquely almost everywhere.

However, due to missingness, the estimating equation in (1) cannot be applied
directly. To formulate the multivariate missingness problem, we further define
the response indicator vector $R=(r_{1},r_{2},\cdots,r_{p})$ for
$Y=(y_{1},y_{2},\cdots,y_{p})$ as

$\displaystyle r_{j}=\left\\{\begin{array}[]{ll}1&\text{if $y_{j}$ is
observed}\\\ 0&\text{otherwise},\end{array}\right.$ (4)

where $j=1,2,\cdots,p$. We assume that the response mechanism is missing at
random in the sense of Rubin (1976). We decompose $Y=(Y_{obs},Y_{mis})$, where
$Y_{obs}$ and $Y_{mis}$ represent the observed and missing parts of $Y$,
respectively. Thus, the missing-at-random assumption is described as

$\displaystyle\mathrm{pr}\left\\{R=(r_{1},r_{2},\cdots,r_{p})\mid
Y_{obs},Y_{mis}\right\\}=\mathrm{pr}\left\\{R=(r_{1},r_{2},\cdots,r_{p})\mid
Y_{obs}\right\\},$ (5)

for any $r_{j}\in\\{0,1\\}$, $j=1,2,\cdots,p$.

Under the missing-at-random assumption, a consistent estimator of $\theta$ can
be obtained by solving the following estimating equation:

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}E\left\\{U(\theta;Y_{i})\mid
Y_{i,obs}\right\\}=0,$ (6)

where it is understood that $E\left\\{U(\theta;Y_{i})\mid
Y_{i,obs}\right\\}=U(\theta;Y_{i})$ if $Y_{i,obs}=Y_{i}$. To compute the
conditional expectation in (6), the parametric fractional imputation (PFI)
method of Kim (2011) can be developed. To apply the parametric fractional
imputation, we can assume that the random vector $Y$ follows a parametric
model in that $F_{0}(Y)\in\left\\{F_{\zeta}(Y):\zeta\in\Omega\right\\}$. In
the parametric fractional imputation, $M$ imputed values for $Y_{i,mis}$, say
$\\{Y_{i,mis}^{*(1)},Y_{i,mis}^{*(2)},\cdots,Y_{i,mis}^{*(M)}\\}$ are
generated from a proposal distribution with the same support of $F_{0}(Y)$ and
are assigned with fractional weights, say
$\left\\{w_{i1}^{*},w_{i2}^{*},\cdots,w_{iM}^{*}\right\\}$, so that a
consistent estimator of $\theta$ can be obtained by solving

$\sum_{i=1}^{n}\sum_{j=1}^{M}w_{ij}^{*}U(\theta;Y_{i,obs},Y_{i,mis}^{*(j)})=0,$

where the fractional weights are constructed to satisfy
$\sum_{j=1}^{M}w_{ij}^{*}U(\theta;Y_{i,obs},Y_{i,mis}^{*(j)})\cong
E\left\\{U\left(\theta;Y_{i}\right)\mid Y_{i,obs}\right\\}$ as closely as
possible, with $\sum_{j=1}^{M}w_{ij}^{*}=1$. In Kim (2011), importance
sampling idea is used to compute the fractional weights.

However, for multivariate missing data, it is not easy to find a joint
distribution family $\left\\{F_{\zeta}(Y):\zeta\in\Omega\right\\}$ correctly.
If the joint distribution family
$\left\\{F_{\zeta}(Y):\zeta\in\Omega\right\\}$ is misspecified, the parametric
fractional imputation can lead to biased inference. All aforementioned
concerns motivate us to consider a more robust fractional imputation method
using Gaussian mixture models, which cover a wider class of parametric models.

##  3 Proposed method

We assume that the random vector $Y$ follows a Gaussian mixture model

$\displaystyle f(Y;\alpha,\zeta)=\sum_{g=1}^{G}\alpha_{g}f(Y;\zeta_{g}),$ (7)

where $G$ is the number of mixture component, $\alpha_{g}\in(0,1)$ is the
mixture proportion satisfying $\sum_{g=1}^{G}\alpha_{g}=1$, and
$f(\cdot;\zeta_{g})$ is the density function of multivariate normal
distribution with parameter $\zeta_{g}=\left\\{\mu_{g},\Sigma_{g}\right\\}$.
Here, we consider the same values of $\Sigma_{g}=\Sigma$ across all components
to get a parsimonious model. The proposed joint model in (7) can be easily
extended to use the group-dependent variance $\Sigma_{g}$, as in Di Zio et al.
(2007).

To formulate the proposal, define the group indicator vector
$Z=(z_{1},z_{2},\cdots,z_{G})$, where $z_{g}=1$ and $z_{j}=0$ for all $j\neq
g$, if sample unit belongs to the $g$-th group. Note that $Z$ is a latent
variable with parameter $\mathrm{pr}(z_{g}=1)=\alpha_{g}$, satisfying
$\sum_{g=1}^{G}\alpha_{g}=1$. We assume that the Gaussian mixture model in (7)
satisfies the strong first-order identifiability assumption (Chen, 1995; Liu
and Shao, 2003; Chen and Khalili, 2008), where the first-order derivatives of
$f(Y;\alpha,\zeta)$ respect to all parameters are linearly independent. Using
$Z$ variable, we can express $f(Y)=\sum_{g=1}^{G}\mathrm{pr}(z_{g}=1)f(Y\mid
z_{g}=1)$, which leads to the marginal distribution in (7).

To handle item nonresponse, we propose to use the fractional imputation method
to impute the missing values. Note that, the joint predictive distribution of
$(Y_{mis},Z)$ given $Y_{obs}$ can be written as $f(Y_{mis},Z\mid
Y_{obs})=f(Z\mid Y_{obs})f(Y_{mis}\mid Y_{obs},Z)$, which implies that the
prediction model for $Y_{mis}$ is

$\displaystyle f(Y_{mis}\mid Y_{obs})=\sum_{g=1}^{G}\mathrm{pr}(z_{g}=1\mid
Y_{obs})f(Y_{mis}\mid Y_{obs},z_{g}=1).$ (8)

The first part in (8) can be obtained by

$\displaystyle\mathrm{pr}(z_{g}=1\mid Y_{obs})=\frac{f(Y_{obs}\mid
z_{g}=1)\alpha_{g}}{\sum_{g=1}^{G}f(Y_{obs}\mid z_{g}=1)\alpha_{g}},$

where $Y_{obs}\mid(z_{g}=1)$ is normal. The second part of (8) , which is
$Y_{mis}\mid(Y_{obs},z_{g}=1)$, is also normal. Therefore, the EM algorithm
for the proposed fractional imputation using Gaussian mixture models (FIGMM)
can be described as follows:

  * I-step: To generate $Y_{i,mis}^{*}$ from $f(Y_{i,mis}\mid Y_{i,obs};\alpha^{(t)},\zeta^{(t)})$ in (8), we use the following two-step method:

    * Step 1: Compute

$\displaystyle p_{ig}^{(t)}=\frac{f(Y_{i,obs}\mid
z_{ig}=1;\zeta_{g}^{(t)})\alpha_{g}^{(t)}}{\sum_{g=1}^{G}f(Y_{i,obs}\mid
z_{ig}=1;\zeta_{g}^{(t)})\alpha_{g}^{(t)}},$

where $f(Y_{i,obs}\mid z_{ig}=1;\zeta_{g})$ is the marginal density of
$Y_{i,obs}$ derived from $(Y_{i,obs},Y_{i,mis})\mid(z_{ig}=1)\sim
N(\mu_{g},\Sigma)$. Generate
$(M_{i1}^{(t)},M_{i2}^{(t)},\cdots,M_{iG}^{(t)})\sim\mathrm{Multinomial}(M;{\mathrm{P}_{i}}^{(t)})$,
where ${\mathrm{P}_{i}}^{(t)}=(p_{i1}^{(t)},\cdots,p_{iG}^{(t)})$ with
$M_{ij}^{(t)}>0$.

    * Step 2: For each $g=1,2,\cdots,G$, we generate $M_{ig}^{(t)}$ independent realizations of $Y_{i,mis}^{*}$, say $\\{Y_{i,mis}^{*(gj)},j=1,2,\cdots,M_{ig}^{(t)}\\}$, from the conditional distribution $f(Y_{i,mis}\mid Y_{i,obs},z_{ig}=1;\zeta_{g}^{(t)})$, which is also normal.

  * W-step: Compute the fractional weights for $Y_{i,mis}^{*(gj)}$ as $w_{igj(t)}^{*}=p_{ig}^{(t)}/M_{ig}^{(t)}$. Note that $\sum_{g=1}^{G}\sum_{j=1}^{M_{ig}^{(t)}}w_{igj(t)}^{*}=1$, for each $i=1,2,\cdots,n$. Using $(w_{igj(t)}^{*},Y_{i,mis}^{*(gj)})$, we can compute

$\displaystyle
Q^{*}(\alpha,\zeta\mid\alpha^{(t)},\zeta^{(t)})=\sum_{i=1}^{n}\sum_{g=1}^{G}\sum_{j=1}^{M_{ig}^{(t)}}w_{igj(t)}^{*}\left\\{\log\alpha_{g}+\log
f(Y_{i}^{*(gj)}\mid z_{ig}=1;\zeta_{g})\right\\},$ (9)

where $Y_{i}^{*(gj)}=(Y_{i,obs},Y_{i,mis}^{*(gj)})$. If $\delta_{i}=1$, then
$Y_{i}^{*(gj)}=Y_{i}$.

  * M-step: Update the parameters by maximizing (9) with respect to $(\alpha,\zeta)$. The solution is

$\displaystyle\alpha_{g}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\frac{M_{ig}^{(t)}}{M};\qquad\mu_{g}^{(t+1)}=\frac{\sum_{i=1}^{n}\sum_{j=1}^{M_{ig}^{(t)}}w_{igj(t)}^{*}Y_{i}^{*(gj)}}{\sum_{i=1}^{n}\sum_{j=1}^{M_{ig}^{(t)}}w_{igj(t)}^{*}};$
$\displaystyle\Sigma^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\sum_{g=1}^{G}\sum_{j=1}^{M_{ig}^{(t)}}w_{igj(t)}^{*}\left(Y_{i}^{*(gj)}-\mu_{g}^{(t+1)}\right)\left(Y_{i}^{*(gj)}-\mu_{g}^{(t+1)}\right)^{\mathrm{\scriptscriptstyle
T}}.$

Repeat I-step to M-step until the convergence is achieved. Then, the final
estimator, say $\hat{\theta}_{FIGMM}$, of $\theta$ can be obtained by solving
the fractionally imputed estimating equation, given by

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\sum_{g=1}^{G}\sum_{j=1}^{M_{ig}}w_{igj}^{*}U(\theta;Y_{i}^{*(gj)})=0,$
(10)

where $w_{igj}^{*}$ are the final fractional weights and $M_{ig}$ are the
final imputation sizes for group $g$, after convergence of the EM algorithm.

###### Remark 1

We now briefly discuss variance estimation of $\hat{\theta}_{FIGMM}$. To
estimate the variance of $\hat{\theta}_{FIGMM}$, replication methods, such as
jackknife, can be used. First note that, the fractional weight assigned to
$Y_{i}^{*(gj)}$ is
$w_{igj}^{*}=\hat{p}_{ig}M_{ig}^{-1}:=\hat{p}_{ig}\hat{\pi}_{2j\mid ig}$,
where $\hat{p}_{ig}$ is obtained from

$\displaystyle\hat{p}_{ig}=\frac{f(Y_{i,obs}\mid
z_{ig}=1;\hat{\zeta}_{g})\hat{\alpha}_{g}}{\sum_{g=1}^{G}f(Y_{i,obs}\mid
z_{ig}=1;\hat{\zeta}_{g})\hat{\alpha}_{g}}.$ (11)

Thus, the $k$-th replicate of $w_{igj}^{*}$ can be obtained by

$\displaystyle w_{igj}^{*(k)}=\hat{p}_{ig}^{(k)}\hat{\pi}_{2j\mid ig}^{(k)},$
(12)

where $\hat{p}_{ig}^{(k)}$ is obtained from (11) using $\hat{\zeta}^{(k)}$ and
$\alpha_{g}^{(k)}$, the $k$-th replicate of $\hat{\zeta}$ and
$\hat{\alpha}_{g}$ respectively, and

$\displaystyle\hat{\pi}_{2j\mid ig}^{(k)}\propto\frac{f(Y_{i,mis}^{*(gj)}\mid
Y_{i,obs},z_{ig}=1;\hat{\zeta}_{g}^{(k)})}{f(Y_{i,mis}^{*(gj)}\mid
Y_{i,obs},z_{ig}=1;\hat{\zeta}_{g})}$

and $\sum_{j=1}^{M_{ig}}\hat{\pi}_{2j\mid ig}^{(k)}=1$. The calculation of
$\hat{\pi}_{2j\mid ig}^{(k)}$ is based on the idea of importance sampling.
Construction of replicate fractional weights using importance sampling idea
has been used in Berg et al. (2016).

The replicate parameter estimates $(\hat{\alpha}^{(k)},\hat{\zeta}^{(k)})$ are
computed by maximizing

$\displaystyle l_{obs}^{(k)}(\alpha,\zeta)=\sum_{i=1}^{n}w_{i}^{(k)}\log
f_{obs}(Y_{i,obs};\alpha,\zeta)$ (13)

respect to $(\alpha,\zeta)$, where
$f_{obs}(Y_{i,obs};\alpha,\zeta)=\sum_{g=1}^{G}\alpha_{g}f(Y_{i,obs}\mid
z_{ig}=1;\zeta_{g})$, and $w_{i}^{(k)}$ is the $k$-th replicate of
$w_{i}=n^{-1}$. The maximizer of $l_{obs}^{(k)}(\alpha,\zeta)$ in (13) can be
obtained by applying the same EM algorithm using replicate weights and
replicate fractional weights in the W-step and M-step. There is no need to
repeat I-step. Variance estimation for $\hat{\theta}_{FIGMM}$ can be obtained
by computing the $k$-th replicate of $\hat{\theta}_{FIGMM}$ from

$\displaystyle\sum_{i=1}^{n}w_{i}^{(k)}\sum_{g=1}^{G}\sum_{j=1}^{M_{ig}}w_{igj}^{*(k)}U(\theta;Y_{i}^{*(gj)})=0.$
(14)

For example, the jackknife variance estimator of $\hat{\theta}_{FIGMM}$ can be
obtained by

$\displaystyle\widehat{\mathrm{var}}\left(\hat{\theta}_{FIGMM}\right)=\frac{n-1}{n}\sum_{i=1}^{n}\left(\hat{\theta}_{FIGMM}^{(k)}-\hat{\theta}_{FIGMM}\right)^{2},$

where $\hat{\theta}_{FIGMM}^{(k)}$ is the $k$-th replicate of
$\hat{\theta}_{FIGMM}$ obtained from (14).

##  4 Asymptotic theory

In our proposed fractional imputation method using Gaussian mixture models in
§3, we have assumed that the size of mixture components, $G$, is known. In
practice , $G$ is often unknown and we need to estimate it from the sample
data. If $G$ is larger than necessary, the proposed mixture model may be
subject to overfitting and increase its variance. If $G$ is small, then the
approximation of the true distribution cannot provide accurate prediction due
to its bias. Hence, we can allow the model complexity parameter $G$ to depend
on the sample size $n$, say $G=G(n)$. The choice of $G$ under complete data
has been well explored in the literature. The popular methods are based on
Bayesian information criterion (BIC) and Akaike’s information criterion (AIC).
See Wallace and Dowe (1999), Windham and Cutler (1992), Schwarz et al. (1978),
Fraley and Raftery (1998), Keribin (2000) and Dasgupta and Raftery (1998). The
alternative way of using SCAD penalty (Fan and Li, 2001) is studied in Chen
and Khalili (2008) and Huang et al. (2017).

In this paper, we consider using the Bayesian information criterion to select
$G$. Under multivariate missingness, we do not have the complete log-
likelihood function. Thus, we use the observed log-likelihood function to
serve the role of the complete log-likelihood function in computing the
information criterion, in the sense that

$\displaystyle\mathrm{BIC}(G)=-2\sum_{i=1}^{n}\log\left\\{\sum_{g=1}^{G}\hat{\alpha}_{g}f(Y_{i,obs}\mid
z_{ig}=1;\hat{\zeta}_{g})\right\\}+(\log n)\phi(G),$ (15)

under the assumption of $\Sigma_{g}=\Sigma$, where
$(\hat{\alpha},\hat{\zeta})$ are the estimators obtained from the proposed
method and $\phi(G)$ is a monotone increasing function of $G$. In (15),
$\phi(G)=G+Gp$ if ignoring constant terms. However, our model selection
framework and theoretical results can be directly applied to any general
penalty function $\phi(G)$. Using Gaussian mixture models, the observed log-
likelihood function is expressed as a closed-form.

In this section, we first establish the consistency of model selection using
(15) under the Gaussian mixture model assumption. After that, we establish
some asymptotic results when the Gaussian mixture model assumption is
violated.

To establish the first part, assume that $\\{Y_{1},Y_{2},\cdots,Y_{n}\\}$ is a
random sample from
$f_{0}(Y)=\sum_{g=1}^{G^{o}}\alpha_{g}^{o}f(Y;\zeta_{g}^{o})$, where
$(G^{o},\alpha^{o},\zeta^{o})$ are true parameter values. For
$\zeta_{g}^{o}=(\mu_{g}^{o},\Sigma^{o})$, we need the following regularity
assumptions:

  * (A1) The mean vectors for each mixture component is bounded uniformly, in the sense of $\|\mu_{g}^{o}\|\leq C_{1}$, for $g=1,2,\cdots,G^{o}$.

  * (A2) $\|\Sigma^{o}\|\leq C_{2}$. Furthermore, $\Sigma^{o}$ is nonsingular.

The first assumption means the first moment is bounded. Assumption (A2) is to
make sure that $\Sigma^{0}$ is bounded and nonsingular. Both assumptions are
commonly used.

To establish the model consistency, we furthermore make the additional
assumptions on the response mechanism:

  * (A3) The response rate for $y_{j}$ is bounded below from 0, say $\lim_{n}n^{-1}\sum_{i=1}^{n}r_{ij}>C_{3}$ , for $j=1,2,\cdots,p$, where $C_{3}>0$ is a constant.

  * (A4) The response mechanism satisfies the mising-at-random condition in (5).

The following theorem shows that the true number of mixture components can be
selected by minimizing $\text{BIC}(G)$ in (15) consistently.

###### Theorem 4.1

Assume the true density $f_{0}$ is the Gaussian mixture model, satisfying
(A1)–(A2). Let $\hat{G}$ be the minimizer of $\mathrm{BIC}(G)$ in (15). Under
assumptions (A3)–(A4), we have

$\displaystyle\mathrm{pr}(\hat{G}=G^{o})\xrightarrow{}1,$

as $n\xrightarrow{}\infty$, where $G^{o}$ is the true number of mixture
components.

The proof of Theorem 4.1 is shown in the Supplementary Material. Theorem 4.1
states that minimizing $\text{BIC}(G)$ consistently selects the true mixture
components under the assumption that the true distribution is in the Gaussian
mixture model.

Now, in the second scenario, the true distribution does not necessary belong
to the class of Gaussian mixture models, Thus, we first establish the
following lemma to measure how well Gaussian mixture model can approximate the
arbitrary density function. We furthermore make additional assumptions about
the true density function $f_{0}$. Use $E_{0}$ to denote the expectation
respect to $f_{0}$.

  * (A5) Assume $f_{0}(Y)$ is continuous with $\int Y^{2}f_{0}(Y)<\infty$.

  * (A6) Assume $E_{0}\left\\{\partial f(Y)/\partial\alpha\right\\}<\infty$ and $E_{0}\left\\{\partial f(Y)/\partial\mu\right\\}<\infty$, where $f(Y)=\sum_{g=1}^{G}\alpha_{g}f(Y;\mu_{g},\Sigma)$. Moreover, assume $E_{0}\left\\{f(Y)^{-2}\right\\}<\infty$.

Assumption (A5) is satisfied for any continuous random variable with bounded
second moments. Assumption (A6) is true for any finite Gaussian mixture model
and $f_{0}$ has a valid moment generating function.

###### Lemma 4.2

Under assumptions (A5)–(A6) and missing at random, for any $\epsilon>0$, there
exist $G=\epsilon^{-\gamma}$, such that

$\displaystyle\|f_{0}-\hat{f}\|_{1}=O(\epsilon),$ (16)
$\displaystyle\mathrm{var}(f_{0}-\hat{f})=O(\epsilon^{-\gamma}n^{-1}),$ (17)

with probability one, where
$\hat{f}(Y)=\sum_{g=1}^{G}\hat{\alpha}_{g}f(Y;\hat{\mu}_{g},\hat{\Sigma})$ is
obtained from the proposed method in §3, $\gamma>0$ and
$\|f_{0}-\hat{f}\|_{1}=\int|f_{0}(Y)-\hat{f}(Y)|f_{0}(Y)\mathrm{d}Y.$

The proof of Lemma 4.2 is presented in the Supplementary Material. If $f_{0}$
is a density function of the Gaussian mixture model, then $\gamma=0$ and by
Theorem 4.1, our proposed $\mathrm{BIC}(G)$ can select the true model
consistently. For any $f_{0}$ satisfies (A5)–(A6) and is not a finite Gaussian
mixture model, the bias can goes to 0 as $G\xrightarrow{}\infty$ from (16).
The variance will increase as $G\xrightarrow{}\infty$ from (17) for fixed $n$.
There is a trade-off between bias and variance for the divergence case
$(\gamma>0,G\xrightarrow{}\infty)$.

Using Lemma 4.2, we can further establish the $\sqrt{n}$-consistency of
$\hat{\theta}_{FIGMM}$. The following assumptions are the sufficient
conditions to obtain the $\sqrt{n}$-consistency.

  * (A7) $E_{0}\left\\{U^{2}(\theta;Y_{i})\right\\}<\infty$.

  * (A8) $\gamma\in(0,2)$.

  * (A9) $\epsilon=O(n^{-1/(2-\Delta)})$, for any $\Delta\in(0,2)$.

###### Theorem 4.3

Under assumptions (A5)–(A9), $\gamma+\Delta<2$ and MAR, we have

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\sum_{g=1}^{G}\sum_{j=1}^{M_{ig}}w_{igj}^{*}U(\theta;Y_{i}^{*(gj)})\cong
J_{1}+o_{p}(n^{-1/2}),$ (18)

where $J_{1}=n^{-1}\sum_{i=1}^{n}E_{0}\left\\{U(\theta;Y_{i})\mid
Y_{i,obs}\right\\}$, if $M=\sum_{i,g}\\{M_{ig}\\}\xrightarrow{}\infty$.
Furthermore, we have

$\displaystyle\sqrt{n}(\hat{\theta}_{FIGMM}-\theta_{0})\xrightarrow{}N(0,\Sigma),$
(19)

for some $\Sigma$ which is positive definite and $\theta_{0}$ satisfies
$E_{0}\left\\{U(\theta_{0};Y)\right\\}=0.$

The proof of (18) is shown in the Supplementary Material and (19) can be
directly derived from (18). From Theorem 4.3, we have
$G=O(n^{\gamma/(2-\Delta)})=o(n)\xrightarrow{}\infty$ with the rate smaller
than $n$. Thus, even under non-Gaussian mixture families, our proposed method
still enjoys $\sqrt{n}$-consistency.

##  5 Extension

In Section 3, we assume that $Y$ is fully continuous. However, in practice,
categorical variables can be used to build imputation models. We extend our
proposed method to incorporate the categorical variable as a covariate in the
model.

To introduce the proposed method, we first introduce the conditional Gaussian
mixture model. Suppose that $(X,Y)$ is a random vector where $X$ is discrete
and $Y$ is continuous. We further assume that $X$ is always observed. To
obtain the conditional Gaussian mixture model, we assume that $Z$ satisfies

$\displaystyle f(Y\mid X,Z)=f(Y\mid Z),$ (20)

in the sense that $Z$ is a partition of the sample such that $Y$ is
homogeneous within each group defined by $Z$. Furthermore, we assume that
$f(Y\mid z_{g}=1)$ follows a Gaussian distribution. Combining these
assumptions, we have the following conditional Gaussian mixture model

$\displaystyle f(Y\mid X)=\sum_{g=1}^{G}\alpha_{g}(X)f(Y\mid z_{g}=1),$ (21)

where $\alpha_{g}(X)=\mathrm{pr}(z_{g}=1\mid X)$ and $f(Y\mid z_{g}=1)$ is the
density function of the normal distribution with parameter
$\zeta_{g}=\left\\{\mu_{g},\Sigma\right\\}$. We also assume that the
identifiability conditions in (21) hold.

Using the argument similar to (8), the predictive model of $Y_{i,mis}$ under
(20) can be expressed as

$\displaystyle f(Y_{i,mis}\mid
Y_{i,obs},X_{i})=\sum_{g=1}^{G}\mathrm{pr}(z_{ig}=1\mid
Y_{i,obs},X_{i})f(Y_{i,mis}\mid Y_{i,obs},z_{ig}=1),$ (22)

where $f(Y_{i,mis}\mid Y_{i,obs},z_{ig}=1)$ can be derived from
$(Y_{i,obs},Y_{i,mis})\mid(z_{ig}=1)\sim N(\mu_{g},\Sigma)$. The posterior
probability of $z_{ig}=1$ given the observed data is

$\displaystyle\mathrm{pr}(z_{ig}=1\mid Y_{i,obs},X_{i})=\frac{f(Y_{i,obs}\mid
z_{ig}=1)\mathrm{pr}(z_{ig}=1\mid X_{i})}{\sum_{g=1}^{G}f(Y_{i,obs}\mid
z_{ig}=1)\mathrm{pr}(z_{ig}=1\mid X_{i})}.$

Therefore, the proposed fractional imputation using conditional Gaussian
mixture models can be summarized as follows:

  * I-step: Creating $M$ imputed values of $Y_{i,mis}$ from (22) can be described in the following two steps.

    * Step 1: For each $g=1,2,\cdots,G$, given the current parameter values $(\alpha_{g}^{(t)},\zeta_{g}^{(t)})$, the posterior probabilities of $z_{ig}=1$ given $(Y_{i,obs},X_{i})$ can be obtained from

$\displaystyle p_{ig}^{(t)}=\frac{f(Y_{i,obs}\mid
z_{ig}=1;\zeta_{g}^{(t)})\alpha_{g}^{(t)}(X_{i})}{\sum_{g=1}^{G}f(Y_{i,obs}\mid
z_{ig}=1;\zeta_{g}^{(t)})\alpha_{g}^{(t)}(X_{i})}.$

    * Step 2: Generate $M$ imputed values of $Y_{i,mis}$ following the same procedure of I-step in §3.

  * W-step: Update the fractional weights for $Y_{i}^{*(gj)}=(Y_{i,obs},Y_{i,mis}^{*(gj)})$ as $w_{igj(t)}^{*}=p_{ig}^{(t)}/{M_{ig}^{(t)}}$, for $j=1,2,\cdots M_{ig}^{(t)}$ and $\sum_{g=1}^{G}M_{ig}^{(t)}=M$. Note that $\sum_{g=1}^{G}\sum_{j=1}^{M_{ig}^{(t)}}w_{igj(t)}^{*}=1$.

  * M-step: Update the parameter values by maximizing

$\displaystyle
Q^{*}(\alpha,\zeta\mid\alpha^{(t)},\zeta^{(t)})=\sum_{i=1}^{n}\sum_{g=1}^{G}\sum_{j=1}^{M_{ig}^{(t)}}w_{igj(t)}^{*}\left\\{\log\alpha_{g}(X_{i})+\log
f(Y_{i}^{*(gj)}\mid z_{ig}=1;\zeta_{g})\right\\},$

respect to $(\alpha,\zeta)$.

Repeat I-step to M-step iteratively until convergence is achieved. The final
estimator of $\theta$ can be obtained by solving the fractionally imputed
estimating equation in (10). Note that the proposed method builds the
proportion vector of mixture components into a function of auxiliary variable
and assumes that the mixture components share the same mean and variance
structure. Thus, the proposed method can borrow information across different
$X$ values. Moreover, the auxiliary information is incorporated to build a
more flexible class of joint distributions.

##  6 Numerical Studies

We consider two simulation studies to evaluate the performance of the proposed
methods. The first simulation study is used to check the performance of the
proposed imputation method using Gaussian mixture models under multivariate
continuous variables. The second simulation study considers the case of
multivariate mixed categorical and continuous variables. To save space, we
only present the first simulation study. The second simulation study is
presented in the Supplementary Material.

In the first simulation study, we consider the following models for generating
$Y_{i}=(Y_{i1},Y_{i2},Y_{i3})$.

  * 1.

M1: A mixture distribution with density
$f(Y)=\sum_{g=1}^{3}\alpha_{g}f_{g}(Y)$, where
$(\alpha_{1},\alpha_{2},\alpha_{3})=(0\mbox{$\cdot$}3,0\mbox{$\cdot$}3,0\mbox{$\cdot$}4)$
and $f_{g}(Y)$ is a density function for multivariate normal distribution with
mean $\mu_{g}$ and variance

$\displaystyle\Sigma(\rho)=\left(\begin{matrix}1&\rho&\rho^{2}\\\
\rho&1&\rho&\\\ \rho^{2}&\rho&1\\\ \end{matrix}\right).$

Let $\rho=0\mbox{$\cdot$}7$ and
$\mu_{1}=(-3,-3,-3,),\mu_{2}=(1,1,1),\mu_{3}=(5,5,5)$.

  * 2.

M2: Use the same model as M2 except for $f_{1}(Y)$, where $f_{1}(Y)$ is a
product of the density for the exponential distribution with rate parameter 1.

  * 3.

M3: $Y_{i1}=1+e_{i1},Y_{i2}=0\mbox{$\cdot$}5Y_{i1}+e_{i2}$ and
$Y_{i3}=Y_{i2}+e_{i3}$, where $e_{i1},e_{i2},e_{i3}$ are independently
generated from $N(0,1)$, $\text{Gamma}(1,1)$ and $\chi_{1}^{2}$ distributions,
respectively.

  * 4.

M4: Generate $(Y_{i1},Y_{i2})$ independently from a Gaussian distribution with
mean $(1,2)$ and variance

$\displaystyle\left(\begin{matrix}1&0\mbox{$\cdot$}5\\\
0\mbox{$\cdot$}5&1\end{matrix}\right).$

Let $Y_{i3}=Y_{i2}^{2}+e_{i3}$, where $e_{i3}\sim N(0,1)$.

In M1, a Gaussian mixture model with $G=3$ is used to generate the samples. A
non-Gaussian mixture distribution is used in M2 to check the robustness of the
imputation methods. M3 and M4 are used to check the performance of the
imputation methods under skewness and nonlinearity, respectively.

The size for each realized sample is $n=500$. Once the complete sample is
obtained, for $y_{ij},j=2,3$, we select 25% of the sample independently to
make missingness with the selection probabilities equal to $\pi_{ij}$, where
$\mathrm{logit}(\pi_{i2})=-0\mbox{$\cdot$}8+0\mbox{$\cdot$}4y_{i1}$,
$\mathrm{logit}(\pi_{i3})=0\mbox{$\cdot$}4-0\mbox{$\cdot$}8y_{i1}$, and
$\mathrm{logit}(u)=\exp(u)/\left\\{1+\exp(u)\right\\}$. Since we assume
$y_{i1}$ are fully observed, the response mechanism is missing at random.

The overall missing rate is approximate 55%. For each realized incomplete
samples, we apply the following methods:

  * [Full]: As a benchmark, we use the full samples to estimate parameters. Confidence intervals with 95% coverage rates are constructed using Wald method.

  * [CC]: Only use the complete cases to estimate parameters and construct confidence intervals.

  * [MICE]: Apply multivariate imputation by chained equations (Buuren and Groothuis-Oudshoorn, 2011). The predictive mean matching is used as a default. The variance estimators are obtained using Rubin’s formula and confidence intervals are constructed using Wald method.

  * [PFI]: Parametric fractional imputation method of Kim (2011), where we assume the joint distribution is a multivariate normal distribution.

  * [SFI]: The proposed semiparametric fractional imputation method using Gaussian mixture models, where the number of components $G$ is selected using the BIC in (15). The confidence interval is implemented using the variance estimator presented in Remark 1.

The parameters of interest are the population means and the population
proportions. For $Y=(y_{1},y_{2},y_{3})$, define
$\theta_{2}=E(y_{2}),\theta_{3}=E(y_{3})$. The parameters of population
proportions are defined as
$P_{2}=\mathrm{pr}(y_{2}<c_{2}),P_{3}=\mathrm{pr}(y_{3}<c_{3})$, where
$(c_{2},c_{3})=(-2,-2)$ for M1 and M2, $(c_{2},c_{3})=(2,3)$ for M3,
$(c_{2},c_{3})=(2,5)$ for M4. The simulation is repeated for $B=2,000$ times.

To evaluate the above methods, the relative mean square error (RMSE) is
defined as

$\displaystyle\mathrm{RMSE}=\frac{\mathrm{MSE}_{\text{method}}}{\mathrm{MSE}_{\text{Full}}}\times
100,$

where $\mathrm{MSE}_{\text{method}}$ is the mean squared error of the
parameters of the current method and $\mathrm{MSE}_{\text{Full}}$ is the mean
squared error of the parameters of using full samples. The simulation results
of RMSE and average coverage probability are presented in Table 1. The
histograms of selected $G$ values using the proposed BIC method are shown in
Figure 1.

Table 1: The relative mean squared errors and coverage probabilities (in parentheses) of the five methods computed from a Monte Carlo simulation of size $B=2,000$. Full, the method uses full samples; CC, the method only uses complete cases; MICE, multiple imputation uses chain equations; PFI, parametric fractional imputation; SFI, the proposed semiparametric fractional imputation method uses Gaussian mixture models. Model | Method | $\theta_{2}$ | $\theta_{3}$ | $P_{2}$ | $P_{3}$
---|---|---|---|---|---
M1 | Full | 100(94$\cdot$9) | 100(95$\cdot$0) | 100(94$\cdot$5) | 100(94$\cdot$9)
CC | 3233(2$\cdot$9) | 3010(3$\cdot$5) | 2639(4$\cdot$5) | 2407(5$\cdot$6)
MICE | 100(95$\cdot$2) | 105(93$\cdot$9) | 102(94$\cdot$7) | 145(90$\cdot$6)
PFI | 100(94$\cdot$6) | 104(93$\cdot$3) | 102(94$\cdot$0) | 141(89$\cdot$1)
SFI | 100(94$\cdot$6) | 106(95$\cdot$0) | 102(95$\cdot$3) | 134(95$\cdot$3)
M2 | Full | 100(94$\cdot$1) | 100(94$\cdot$6) | 100(94$\cdot$1) | 100(94$\cdot$5)
CC | 2841(3$\cdot$0) | 2863(3$\cdot$6) | 2284(4$\cdot$8) | 2233(5$\cdot$7)
MICE | 102(94$\cdot$4) | 109(93$\cdot$7) | 103(93$\cdot$8) | 189(86$\cdot$3)
PFI | 102(93$\cdot$8) | 108(94$\cdot$1) | 102(93$\cdot$3) | 180(83$\cdot$8)
SFI | 106(95$\cdot$1) | 107(95$\cdot$4) | 103(93$\cdot$8) | 177(92$\cdot$8)
M3 | Full | 100(94$\cdot$7) | 100(94$\cdot$7) | 100(95$\cdot$2) | 100(95$\cdot$1)
CC | 157(77$\cdot$9) | 155(85$\cdot$6) | 145(77$\cdot$0) | 145(86$\cdot$6)
MICE | 100(94$\cdot$8) | 100(93$\cdot$5) | 117(92$\cdot$5) | 128(90$\cdot$4)
PFI | 100(94$\cdot$6) | 100(95$\cdot$2) | 117(89$\cdot$9) | 127(89$\cdot$7)
SFI | 97(95$\cdot$1) | 83(92$\cdot$8) | 106(93$\cdot$8) | 91(94$\cdot$3)
M4 | Full | 100(94$\cdot$9) | 100(95$\cdot$2) | 100(94$\cdot$3) | 100(95$\cdot$2)
CC | 386(81$\cdot$8) | 353(86$\cdot$2) | 317(86$\cdot$4) | 308(86$\cdot$6)
MICE | 110(94$\cdot$8) | 114(95$\cdot$2) | 128(93$\cdot$5) | 207(86$\cdot$7)
PFI | 108(94$\cdot$8) | 111(94$\cdot$2) | 129(91$\cdot$6) | 197(84$\cdot$4)
SFI | 126(95$\cdot$0) | 124(94$\cdot$7) | 135(94$\cdot$5) | 139(94$\cdot$2)
Figure 1: The histograms of selected $G$ values using the proposed BIC.

Table 1 presents the relative mean squared errors and the corresponding
average coverage probabilities from the above simulation study of size
$B=2,000$. Under M1, the joint model is a Gaussian mixture. The proposed
method obtains almost the same RMSEs for estimating $\theta_{2},\theta_{3}$
and $P_{2}$ with MICE and PFI, but outperforms MICE and PFI for estimating
$P_{3}$. Moreover, coverage probabilities of MICE and PFI are less than 95% in
estimating $\theta_{3}$ and $P_{3}$. SFI uniformly achieves approximate 95%
coverage probabilities. Thus, we can conclude that MICE and PFI are biased,
due to model misspecification. The histogram in Figure 1 shows that most of
selected G values are 3, which is the true number of mixture components.

In M2, instead of Gaussian mixture models, one component is the exponential
distribution. Table 1 shows that SFI outperforms MICE and PFI in term of RMSEs
for estimating proportions and obtains similar performance for estimating
means. The coverage rates for MICE and PFI for $P_{3}$ are poor.

The joint distribution M3 is a skewed distribution. From Tables 1, we can see
that SFI outperforms MICE and PFI uniformly. Furthermore, SFI provides better
coverage probabilities than MICE and PFI for $P_{2}$ and $P_{3}$.

The joint distribution in M4 has a nonlinear mean structure of $Y_{i3}$. Under
M4, SFI obtains the much smaller RMSE than MICE and PFI in estimating $P_{4}$,
but larger RMSE in estimating $\theta_{2},\theta_{3},P_{3}$. However, SFI
achieves consistent confidence intervals with approximate 95% coverage
probabilities. MICE and PFI are biased in interval estimation and coverage
probabilities are much less than 95%. Overall, the performance of SFI is much
better than MICE or PFI in terms of coverage probabilities.

Interestingly, the imputed estimators are sometimes more efficient than the
full sample estimators. This phenomenon, called superefficiency (Meng, 1994),
can happen， when the method-of-moment is used in the full sample estimator.
Yang and Kim (2016) give a rigorous theoretical justification for this
phenomenon.

##  7 Application

In this section, we apply the proposed method in §3 to a synthetic data that
mimics monthly retail trade survey data at U.S. Census Bureau. The synthetic
monthly retail trade survey data was made for contest in a conference and can
be found in http://www.portal-stat.admin.ch/ices5/imputation-contest/. The
sampling scheme is a stratified simple random sample without replacement
sample with six strata: one certain (take-all) and five non-certainty strata.
The sample sizes are computed using Neyman allocation. An overview of the
monthly retail trade survey data is presented in Figure 2.

Figure 2: Overview for log-values of the monthly retail trade survey
data:“mos” is frame measure of size; “Sales00” denotes current month sales for
unit (subject to missing); “asales00” is current month administrative data
value for sales; “Sales01” means prior month sales for unit; “Inventories00”
is current month inventories for unit (Subject to missing); “ainventories00”
is current month administrative data value for inventories; “Inventories01” is
prior month inventories for unit.

The overall response rate is approximately 71%. Current month sales and
inventories are subject to missingness. From Figure 2, we can find that this
monthly retail trade data are highly skewed. From the normal quantile-quantile
plot, normality assumption is violated and there exist three extreme outliers.

To impute current month sales and inventories, we applied the proposed
fractional imputation method and MICE. After implementation, MICE failed to
converge due to high correlations among the survey items. Therefore, we only
present the final analysis results using the proposed method. The final
results are shown in Table 2.

Table 2: Imputation results for the monthly retail trade survey data. Parameter estimation, 95% confidence interval and true values are presented. Estimate, the estimators provided by the proposed method; Truth, the true parameter values from U.S. Census Bureau. Parameter | Estimate | 95% Confidence interval | Truth
---|---|---|---
Mean of Sales00 $(\times 10^{-6})$ | 2$\cdot$28 | (2$\cdot$10, 2$\cdot$46) | 2$\cdot$30
Mean of Inventories00 $(\times 10^{-6})$ | 4$\cdot$76 | (4$\cdot$42, 5$\cdot$10) | 4$\cdot$81
Correlation of Sales00 and Inventories00 | 0$\cdot$97 | (0$\cdot$94, 0$\cdot$99) | 0$\cdot$97

Comparing with the true population statistics, provided by U.S. Census Bureau,
we can see that our proposed fractional imputation method with Gaussian
mixture models works well to preserve the correlation structure and handle
skewness and outliers. In Table 2, we can see that all 95% confidence
intervals contain their true values.

##  8 Discussion

Fractional imputation has been proposed as a tool for frequentist imputation,
as an alternative to multiple imputation. Multiple imputation using Rubin’s
formula can be biased when the model is uncongenial or the point estimator is
not self-efficient (Meng, 1994; Yang and Kim, 2016). In this paper, we have
proposed a semiparametric fractional imputation method using Gaussian mixture
models to handle arbitrary multivariate missing data. The proposed method
automatically selects the size of mixture components and provides a unified
framework for robust imputation. Even if the group size $G$ increases with the
sample size $n$, the resulting estimator enjoys $\sqrt{n}$-consistency. We
have also extended the proposed method to incorporate categorical auxiliary
variable. The flexible model assumption and efficient computation are the main
advantages of our proposed method. An extension of the proposed method to
survey data is a topic of future research. An R software package for the
proposed method is under development.

## References

  * Bacharoglou (2010) Bacharoglou, A. (2010).  Approximation of probability distributions by convex mixtures of Gaussian measures.  Proceedings of the American Mathematical Society 138(7), 2619–2628.
  * Bartlett et al. (2015) Bartlett, J. W., S. R. Seaman, I. R. White, J. R. Carpenter, and A. D. N. Initiative (2015).  Multiple imputation of covariates by fully conditional specification: accommodating the substantive model.  Statistical Methods in Medical Research 24(4), 462–487.
  * Berg et al. (2016) Berg, E., J.-K. Kim, and C. Skinner (2016).  Imputation under informative sampling.  Journal of Survey Statistics and Methodology 4(4), 436–462.
  * Buuren and Groothuis-Oudshoorn (2011) Buuren, S. and K. Groothuis-Oudshoorn (2011).  MICE: Multivariate imputation by chained equations in R.  Journal of Statistical Software 45(3), 1–67.
  * Chen (2010) Chen, H. Y. (2010).  Compatibility of conditionally specified models.  Statistics & Probability Letters 80(7-8), 670–677.
  * Chen (1995) Chen, J. (1995).  Optimal rate of convergence for finite mixture models.  The Annals of Statistics, 221–233.
  * Chen and Khalili (2008) Chen, J. and A. Khalili (2008).  Order selection in finite mixture models with a nonsmooth penalty.  Journal of the American Statistical Association 103(484), 1674–1683.
  * Cheng (1994) Cheng, P. E. (1994).  Nonparametric estimation of mean functionals with data missing at random.  Journal of the American Statistical Association 89(425), 81–87.
  * Dasgupta and Raftery (1998) Dasgupta, A. and A. E. Raftery (1998).  Detecting features in spatial point processes with clutter via model-based clustering.  Journal of the American Statistical Association. 93(441), 294–302.
  * Di Zio et al. (2007) Di Zio, M., U. Guarnera, and O. Luzi (2007).  Imputation through finite gaussian mixture models.  Computational Statistics & Data Analysis 51(11), 5305–5316.
  * Elliott and Stettler (2007) Elliott, M. R. and N. Stettler (2007).  Using a mixture model for multiple imputation in the presence of outliers: the ‘healthy for life’project.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 56(1), 63–78.
  * Fan and Li (2001) Fan, J. and R. Li (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association. 96(456), 1348–1360.
  * Fraley and Raftery (1998) Fraley, C. and A. E. Raftery (1998).  How many clusters? which clustering method? answers via model-based cluster analysis.  The di2007imputationurnal 41(8), 578–588.
  * Huang et al. (2017) Huang, T., H. Peng, and K. Zhang (2017).  Model selection for Gaussian mixture models.  Statistica Sinica 27(1), 147–169.
  * Judkins et al. (2007) Judkins, D., T. Krenzke, A. Piesse, Z. Fan, and W.-C. Haung (2007).  Preservation of skip patterns and covariate structure through semi-parametric whole questionnaire imputation.  In Proceedings of the Section on Survey Research Methods of the American Statistical Association, pp. 3211–3218.
  * Keribin (2000) Keribin, C. (2000).  Consistent estimation of the order of mixture models.  Sankhya A, 49–66.
  * Kim et al. (2014) Kim, H. J., J. P. Reiter, Q. Wang, L. H. Cox, and A. F. Karr (2014).  Multiple imputation of missing or faulty values under linear constraints.  Journal of Business & Economic Statistics 32(3), 375–386.
  * Kim (2011) Kim, J. K. (2011).  Parametric fractional imputation for missing data analysis.  Biometrika 98(1), 119–132.
  * Kim and Fuller (2004) Kim, J. K. and W. Fuller (2004).  Fractional hot deck imputation.  Biometrika 91(3), 559–578.
  * Kim and Shao (2013) Kim, J. K. and J. Shao (2013).  Statistical Methods for Handling Incomplete Data.  CRC Press.
  * Little and Rubin (2002) Little, R. J. and D. B. Rubin (2002).  Statistical Analysis with Missing Data.  John Wiley & Sons.
  * Liu et al. (2013) Liu, J., A. Gelman, J. Hill, Y.-S. Su, and J. Kropko (2013).  On the stationary distribution of iterative imputations.  Biometrika 101(1), 155–173.
  * Liu and Shao (2003) Liu, X. and Y. Shao (2003).  Asymptotics for likelihood ratio tests under loss of identifiability.  The Annals of Statistics 31(3), 807–832.
  * McLachlan and Peel (2004) McLachlan, G. and D. Peel (2004).  Finite Mixture Models.  John Wiley & Sons.
  * Meng (1994) Meng, X.-L. (1994).  Multiple-imputation inferences with uncongenial sources of input.  Statistical Science 9(4), 538–558.
  * Murray and Reiter (2016) Murray, J. S. and J. P. Reiter (2016).  Multiple imputation of missing categorical and continuous values via Bayesian mixture models with local dependence.  Journal of the American Statistical Association 111(516), 1466–1479.
  * Raghunathan et al. (2001) Raghunathan, T. E., J. M. Lepkowski, J. Van Hoewyk, and P. Solenberger (2001).  A multivariate technique for multiply imputing missing values using a sequence of regression models.  Survey Methodology 27(1), 85–96.
  * Rubin (1976) Rubin, D. B. (1976).  Inference and missing data.  Biometrika 63(3), 581–592.
  * Rubin (1996) Rubin, D. B. (1996).  Multiple imputation after 18+ years.  Journal of the American Statistical Association. 91(434), 473–489.
  * Schwarz et al. (1978) Schwarz, G. et al. (1978).  Estimating the dimension of a model.  The Annals of Statistics 6(2), 461–464.
  * Tanner and Wong (1987) Tanner, M. A. and W. H. Wong (1987).  The calculation of posterior distributions by data augmentation.  Journal of the American Statistical Association. 82(398), 528–540.
  * Wallace and Dowe (1999) Wallace, C. S. and D. L. Dowe (1999).  Minimum message length and kolmogorov complexity.  The Computer Journal 42(4), 270–283.
  * Wang and Chen (2009) Wang, D. and S. X. Chen (2009).  Empirical likelihood for estimating equations with missing values.  The Annals of Statistics 37(1), 490–517.
  * Windham and Cutler (1992) Windham, M. P. and A. Cutler (1992).  Information ratios for validating mixture analyses.  Journal of the American Statistical Association. 87(420), 1188–1192.
  * Yang and Kim (2016) Yang, S. and J. K. Kim (2016).  A note on multiple imputation for method of moments estimation.  Biometrika 103(1), 244–251.