# Matrix Completion for Survey Data Prediction with Multivariate Missingness

Xiaojun Mao, Zhonglei Wang and Shu Yang School of Data Science, Fudan
University, Shanghai 200433, P.R.C. Email<EMAIL_ADDRESS>Key
Laboratory of Econometrics, Wang Yanan Institute for Studies in Economics and
School of Economics, Xiamen University, Xiamen, Fujian 361005, P.R.C. Email:
<EMAIL_ADDRESS>of Statistics, North Carolina State University,
North Carolina 27695, U.S.A. Email<EMAIL_ADDRESS>
###### Abstract

The National Health and Nutrition Examination Survey (NHANES) studies the
nutritional and health status over the whole U.S. population with
comprehensive physical examinations and questionnaires. However, survey data
analyses become challenging due to inevitable missingness in almost all
variables. In this paper, we develop a new imputation method to deal with
multivariate missingness at random using matrix completion. In contrast to
existing imputation schemes either conducting row-wise or column-wise
imputation, we treat the data matrix as a whole which allows exploiting both
row and column patterns to impute the missing values in the whole data matrix
at one time. We adopt a column-space-decomposition model for the population
data matrix with easy-to-obtain demographic data as covariates and a low-rank
structured residual matrix. A unique challenge arises due to lack of
identification of parameters in the sample data matrix. We propose a
projection strategy to uniquely identify the parameters and corresponding
penalized estimators, which are computationally efficient and possess desired
statistical properties. The simulation study shows that the doubly robust
estimator using the proposed matrix completion for imputation has smaller mean
squared error than other competitors. To demonstrate practical relevance, we
apply the proposed method to the 2015-2016 NHANES Questionnaire Data.

Key Words: Double robustness; Imputation; Low rank matrix; Missingness at
random.

##  1 Introduction

Survey data are the gold-standard for estimating finite population parameters
and providing a comprehensive overview of the finite population at a given
time. The National Health and Nutrition Examination Survey (NHANES,
https://www.cdc.gov/nchs/nhanes), for example, is a program of studies to
assess the health and nutrition status of the adults and children in the
United States. The survey combines physical examinations and questionnaires
and therefore can be used to provide a thorough and detailed health status
assessment. In the 2015-2016 Questionnaire Data, there are about 39 blocks of
questions, such as dietary behavior and alcohol use, and each block contains
about ten relevant questions.

However, survey data analyses become challenging due to inevitable
multivariate missingness, leading to complex “swiss cheese” patterns. This
occurs due to item nonresponse, when individuals provide answers to partial
but not all questions. Moreover, the missingness rates vary across questions
and are extremely low for sensitive questions such as income. In the the
NHANES 2015-2016 Questionnaire Data, the average and standard error of the
missingness rates are about 0.62 and 0.38, respectively. This phenomenon is
not an exception but a rule for large surveys in the United States, including
the American National Election Studies, American Housing Survey and Current
Population Survey. Inference ignoring the nonresponse items may be
questionable (Rubin, 1976)

Imputation is widely used to handle item nonresponse, and existing methods for
multivariate missingness can be categorized into two types: row-wise
imputation and column-wise imputation. For example, multiple imputation
(Rubin, 1976; Clogg et al., 1991; Fay, 1992; Meng, 1994; Wang and Robins,
1998; Nielsen, 2003; Kim et al., 2006; Kim, 2011; Yang and Kim, 2016) can be
viewed as a row-wise imputation method, which models the joint distribution of
all variables and generates the imputations based on a posterior predictive
distribution of the nonresponse items given the observed ones. However,
multiple imputation is sensitive to model misspecification, especially when
there are a lot of questions subject to non-response. Moreover, it is
computationally intensive, and it quickly becomes infeasible to implement as
the number of questions subject to missingness increases. On the other hand,
hot deck imputation (Chen and Shao, 2000; Kim and Fuller, 2004; Fuller and
Kim, 2005; Andridge and Little, 2010) can be viewed as a column-wise
imputation method. For subject $i$ with missing $y_{ij}$ of the $j$th
question, hot deck imputation methods search among the units with responses to
the $j$th question (referred to as donors for the $j$th question), and impute
the missing $y_{ij}$ by the response from its nearest neighbor based on a
certain distance metric.

In contrast to most existing methods that use either models or distance, we
treat the data matrix as a whole and propose using matrix completion (Candès
and Recht, 2009; Keshavan, Montanari and Oh, 2009; Mazumder, Hastie and
Tibshirani, 2010; Koltchinskii, Lounici and Tsybakov, 2011; Negahban and
Wainwright, 2012; Cai and Zhou, 2016; Robin et al., 2019) as a tool for
imputaiton, which allows exploiting both row and column patterns to impute the
missing values in the whole data matrix at one time. Because there exist
variables that are fully observed, we adopt a column-space-decomposition model
(Mao, Chen and Wong, 2019) for the population data matrix with easy-to-
obtained demographic data as covariates and a low-rank structured residual
matrix. The low-rank structure is due to underlying clusters of individuals
and blocks of questions (Candès and Recht, 2009; van der Linden and Hambleton,
2013; Davenport and Romberg, 2016; Robin et al., 2019). Most works in the
matrix completion literature assume uniform missingness (or equivalently
missingness completely at random), which however is unlikely to hold in the
survey data context. Following Mao, Wong and Chen (2019), we assume that the
missing data mechanism is missingness at random (MAR; Rubin, 1976). Even
though the population risk function identifies the parameter and data matrix
uniquely, the sample risk function lacks identification in general. We propose
a projection strategy so that the new set of parameters can be identified
after projection based on the sample data. For estimation, we consider a risk
function weighted by both design weights and inverse of the estimated response
probabilities, with the nuclear norm to encourage the low-rankness and two
Frobenious norms to improve numerical performance of the penalized estimators.
After imputation for the sample data, we use a doubly robust estimator for the
population means (Kott, 1994; Bang and Robins, 2005; Kim and Haziza, 2014;
Haziza and Rao, 2006; Kang and Schafer, 2007; Kott and Chang, 2010), which is
unbiased when the response model is correctly specified.

The proposed method achieves the following advantages. First, it is
computationally easy. Based on the column-space-decomposition model, we have
modified the objective function so that we can obtain a closed-form solution
to recover the sample data matrix, and only one singular value decomposition
(SVD) of an $n\times L$ matrix is required for computation. Second, it is a
multi-purpose imputation method; that is, a single-imputation system can be
applied to all the survey questions. This is particularly attractive for a
comprehensive analysis of the whole survey data. Third, comparing to fully
parametric methods, we only require a low-rank assumption without any further
specification. For theoretical investigation, we provide regularity conditions
and the asymptotic bounds of the penalized estimators and the doubly robust
estimator.

The paper proceeds as follows. Section 2 provides the basic setup and
estimation procedure of the proposed method. Section 3 discusses the
theoretical properties of the proposed method. A simulation study is conducted
in Section 4 to illustrate the advantage of the proposed method compared with
other competitors. Section 5 presents an application to the NHANES 2015-2016
Questionnaire Data. Some concluding remarks are given in Section 6.

##  2 Basic Setup

###  2.1 Notation, Assumption and Model

Consider a finite population of $N$ subjects with study variables
$U_{N}=\\{(\bm{x}_{i},\bm{y}_{i}):i=1,\ldots,N\\}$. Organize the finite
population in matrix forms $\mathbf{X}_{N}=(x_{ij})\in\mathbb{R}^{N\times d}$
and $\mathbf{Y}_{N}=(y_{ij})\in\mathbb{R}^{N\times L}$, where $\mathbf{X}_{N}$
is fully observed, and $\mathbf{Y}_{N}$ is subject to missingness. We are
interested in estimating $\theta_{j}=N^{-1}\sum_{i=1}^{N}y_{ij}$ for
$j=1,\ldots,L$.

Assume that the finite population is a realization of an infinite population,
called a super-population, and consider a super-population model $\zeta$,

$\mathbf{Y}_{N}=\mathbf{A}_{N}+\bm{\epsilon}_{N},$ (1)

where $\mathbf{A}_{N}\in\mathbb{R}^{N\times L}$ represents the structural
component of the data matrix, and
$\bm{\epsilon}_{N}=(\epsilon_{ij})\in\mathbb{R}^{N\times L}$ is a matrix of
independent errors with $E(\epsilon_{ij})=0$ for $i=1,\ldots,N$ and
$j=1,\ldots,L$. Following Candès and Recht (2009) ,van der Linden and
Hambleton (2013), Davenport and Romberg (2016) and Robin et al. (2019), we
assume that $\mathbf{A}_{N}$ has a low-rank structure, which is reasonable in
the survey context. On the one hand, the finite population can be divided into
groups by demographics such as age, gender, address and occupations. On the
other hand, survey questions can also be grouped into several blocks. For
example, in the NHANES 2015-1016 Questionnaire Data, there exist different
blocks of questions, such as health and nutrition status, education, income
level and so on, and each block contains several relevant questions.

To further incorporate $\mathbf{X}_{N}$ into the model (1), following Mao,
Chen and Wong (2019), we adopt the column-space-decomposition model,

$\mathbf{A}_{N}=\mathbf{X}_{N}\bm{\beta}^{\ast}+\mathbf{B}_{N}^{\ast},$ (2)

where $\bm{\beta}^{\ast}=(\beta_{ij})$ is a $d\times L$ coefficient matrix,
and $\mathbf{B}_{N}^{\ast}$ is an $N\times L$ low-rank matrix, inherited from
the low-rank structure of $\mathbf{A}_{N}$. To avoid identification issues, we
assume that $\mathbf{X}_{N}^{\rm T}\mathbf{B}_{N}^{\ast}=\mathbf{0}$. We also
assume that the elements in the matrices are indexed by $N$ implicitly.

Following Fay (1991) and Shao and Steel (1999), we first have a census with
nonrespondents. Denote $\mathbf{R}_{N}=(r_{ij})\in\mathbb{R}^{N\times L}$ as
the response indicator matrix with $r_{ij}=1$ if $y_{ij}$ is observed and
$r_{ij}=0$ otherwise. Following most of the missing data literature, we assume
that the missing data mechanism is MAR. Specifically, assume that
$\mathbf{X}_{N}$ explains the missing mechanism well in the sense that the
values in $\bm{y}_{i}$ are MAR conditional on $\bm{x}_{i}$. Under MAR, the
response probability becomes
$p_{ij}=\Pr(r_{ij}=1\mid\mathbf{X}_{N},\mathbf{Y}_{N})=\Pr(r_{ij}=1\mid\bm{x}_{i})$.
For regularity reasons, we require $p_{ij}$ to be bounded away from 0 and 1
for all $i$ and $j$. Following most of the empirical literature, we assume
that the response probability follows a logistic regression model,

$p_{ij}=p_{ij}\left(\bm{x}_{i}\right)=\frac{\exp\left\\{\left(1,\bm{x}_{i}^{\rm
T}\right)\bm{\gamma}_{.j}\right\\}}{1+\exp\left\\{\left(1,\bm{x}_{i}^{\rm
T}\right)\bm{\gamma}_{.j}\right\\}},$ (3)

where $\bm{\gamma}_{.j}\in\mathbb{R}^{d+1}$ is the parameter vector specific
for the $j$th column of $\mathbf{Y}_{N}$. Denote
$\mathbf{P}_{N}^{\dagger}=(p_{ij}^{-1})\in\mathbb{R}^{N\times L}$ as the
matrix of the inverse response probabilities.

To motivate the proposed method, we first consider the population data. Denote
$\mathcal{C}(\mathbf{X})$ as the column space of a matrix $\mathbf{X}$ and
$\mathcal{N}(\mathbf{X})=\\{\mathbf{M}\in\mathbb{R}^{N\times
L}:\mathbf{X}^{\rm T}\mathbf{M}=\mathbf{0}\\}$. Under Model (2) and MAR, for
any $\bm{\beta}\in\mathbb{R}^{d\times L}$ and
$\mathbf{B}\in\mathcal{N}(\mathbf{X}_{N})$, the population risk function
$R(\bm{\beta},\mathbf{B})$ is

$R(\bm{\beta},\mathbf{B})=\frac{1}{NL}E\left\|\mathbf{R}_{N}\circ\mathbf{P}_{N}^{\dagger}\circ\mathbf{Y}_{N}-\mathbf{X}_{N}\bm{\beta}-\mathbf{B}\right\|_{F}^{2},$
(4)

where “$\circ$” is the Hardamard product, and
$\|\mathbf{M}\|_{F}=(\sum_{i=1}^{N}\sum_{j=1}^{L}m_{ij}^{2})^{1/2}$ is the
Frobenius norm of an $N\times L$ matrix $\mathbf{M}=(m_{ij})$. Then,
$(\bm{\beta}^{\ast},\mathbf{B}_{N}^{\ast})$ in (2) uniquely minimizes the
population risk function $R(\bm{\beta},\mathbf{B})$; see Mao, Chen and Wong
(2019) for details.

In practice, it is both time-consuming and expensive to conduct a census for a
finite population. Survey sampling has been the gold standard to estimate
finite population parameters based on a relatively small probability sample.
Assume that a sample of size $n$ is selected by a probability sampling design
(Fuller, 2009). Denote $I_{i}$ as the sampling indicator; specifically,
$I_{i}=1$ if the $i$th subject is sampled and 0 otherwise. Let
$\pi_{i}=E(I_{i}\mid U_{N})$ be the inclusion probability of the $i$th
subject, where the expectation is taken with respect to the probability
sampling mechanism. For example, Poisson sampling generates a sample using $N$
independent Bernoulli trials, where $I_{i}$ is generated from a Bernoulli
distribution with success probability $\pi_{i}$ for $i=1,\ldots,N$. Without
loss of generality, assume that the first $n$ subjects of the finite
population are sampled. In the following, without ambiguity, we use
$\mathbf{M}_{n}$ to denote the sample data matrix of the first $n$ rows of a
population data matrix $\mathbf{M}_{N}$. If $\mathbf{Y}_{n}$ is fully
observed, we can obtain a Horvitz-Thompson estimator (Horvitz and Thompson,
1952) of $\theta_{j}$,

$\widehat{\theta}_{j}=\frac{1}{N}\sum_{i=1}^{N}\frac{I_{i}}{\pi_{i}}y_{ij}.$
(5)

It follows that $\widehat{\theta}_{j}$ is a design-unbiased estimator of
$\theta_{j}$, that is, $E(\widehat{\theta}_{j}\mid U_{N})=\theta_{j}$.

In the presence of missingness in $\mathbf{Y}_{n}$, we propose to use matrix
completion as an imputation method to recover the missing values. With the
sample data matrices, we use the following empirical risk function to
approximate the population risk $R(\bm{\beta},\mathbf{B})$ in (4):

$\displaystyle R^{\ast}(\bm{\beta},\mathbf{B}_{n})$ $\displaystyle=$
$\displaystyle\frac{1}{NL}\sum_{i=1}^{N}\frac{I_{i}}{\pi_{i}}\sum_{j=1}^{L}\left\\{\frac{r_{ij}}{p_{ij}}y_{ij}-(\mathbf{X}_{N}\bm{\beta})_{ij}-b_{ij}\right\\}^{2}$
(6) $\displaystyle=$
$\displaystyle\frac{1}{NL}\left\|\mathbf{D}_{n}^{-1/2}\left(\mathbf{R}_{n}\circ\mathbf{P}_{n}^{\dagger}\circ\mathbf{Y}_{n}-\mathbf{X}_{n}\bm{\beta}-\mathbf{B}_{n}\right)\right\|_{F}^{2},$
(7)

where $\mathbf{D}_{n}=\mbox{diag}(\pi_{1},\ldots,\pi_{n})$ is a diagonal
matrix with $\pi_{i}$ being the $(i,i)$th entry. If the sampling mechanism is
non-informative, there is no need to adjust sampling weights for estimating
$\bm{\beta}^{\ast}$ and $\mathbf{B}_{N}^{\ast}$ in (2). Adjusting for sampling
weights, however, achieves two goals. First, the expectation of (6) is the
population risk function $R(\bm{\beta},\mathbf{B})$, so we target for
estimating the population data matrix instead of the sample data matrix.
Second, it allows for informative sampling. Under informative sampling, the
empirical risk function without sampling weights is biased of the population
risk function.

###  2.2 Non-identifiability of $(\bm{\beta}^{\ast},\mathbf{B}_{n}^{\ast})$

In the population risk function (4), $\mathbf{X}_{N}^{\rm
T}\mathbf{B}_{N}^{\ast}=\mathbf{0}$ guarantees that
$(\bm{\beta}^{\ast},\mathbf{B}_{N}^{\ast})$ is identifiable; see Mao, Chen and
Wong (2019). Moreover, the decomposition of $\mathbf{A}_{N}$ into
$\mathbf{X}_{N}\bm{\beta}^{\ast}\in\mathcal{C}(\mathbf{X}_{N})$ and
$\mathbf{B}_{N}^{\ast}\in\mathcal{N}(\mathbf{X}_{N})$ gives benefits for
showing theoretical properties of the estimators and encourages an efficient
algorithm allowing for a closed-form solution of $\widehat{\mathbf{B}}_{N}$.
However, the same decomposition technique may fail to guarantee identification
of parameters in the sample risk function
$R^{\ast}(\bm{\beta},\mathbf{B}_{n})$ in (7) because
$(\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n})^{\rm
T}(\mathbf{D}_{n}^{-1/2}\mathbf{B}_{n})=\mathbf{X}_{n}^{\rm
T}\mathbf{D}_{n}^{-1}\mathbf{B}_{n}$ may not be a zero matrix. Even for simple
random sampling with $\pi_{i}=n/N$ for $i=1,\ldots,N$, we cannot ensure
$\mathbf{X}_{n}^{\rm
T}\mathbf{D}_{n}^{-1}\mathbf{B}_{n}=Nn^{-1}\mathbf{X}_{n}^{\rm
T}\mathbf{B}_{n}=\mathbf{0}$. It means that there is no space restriction for
both $\bm{\beta}$ and $\mathbf{B}_{n}$ in
$R^{\ast}(\bm{\beta},\mathbf{B}_{n})$. Thus, for any
$(\bm{\beta},\mathbf{B}_{n})$ and nonzero $\bm{\beta}_{1}$, we always have
$R^{\ast}(\bm{\beta},\mathbf{B}_{n})=R^{\ast}(\bm{\beta}+\bm{\beta}_{1},\mathbf{B}_{n}-\mathbf{X}_{n}\bm{\beta}_{1})$.

To deal with the lack of identifiability, we consider a decomposition of
$\mathbf{D}_{n}^{-1/2}(\mathbf{R}_{n}\circ\mathbf{P}_{n}^{\dagger}\circ\mathbf{Y}_{n}-\mathbf{X}_{n}\bm{\beta}-\mathbf{B}_{n})$
by

$\mathbf{D}_{n}^{-1/2}\left(\mathbf{R}_{n}\circ\mathbf{P}_{n}^{\dagger}\circ\mathbf{Y}_{n}\right)-\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}\bm{\beta}-\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}(\mathbf{D}_{n}^{-1/2}\mathbf{B}_{n})-\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}^{\perp}(\mathbf{D}_{n}^{-1/2}\mathbf{B}_{n}).$

where
$\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}=\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}(\mathbf{X}_{n}^{\rm
T}\mathbf{D}_{n}^{-1}\mathbf{X}_{n})^{-1}\mathbf{X}_{n}^{\rm
T}\mathbf{D}_{n}^{-1/2}$,
$\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}^{\perp}=\mathbf{I}-\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}$
and $\mathbf{I}$ is the $n\times n$ identity matrix. Denote

$\bm{\beta}^{\ast\prime}=\bm{\beta}^{\ast}+(\mathbf{X}_{n}^{\rm
T}\mathbf{D}_{n}^{-1}\mathbf{X}_{n})^{-1}\mathbf{X}_{n}^{\rm
T}\mathbf{D}_{n}^{-1}\mathbf{B}_{n}^{\ast}\qquad\text{and}\qquad\mathbf{B}_{n}^{\ast\prime}=\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}^{\perp}(\mathbf{D}_{n}^{-1/2}\mathbf{B}_{n}^{\ast}),$

respectively. Then, we have
$\mathbf{B}_{n}^{\ast\prime}\in\mathcal{N}(\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n})$,
so we can decompose the objective function
$R^{\ast}(\bm{\beta},\mathbf{B}_{n})$ as

$\displaystyle
R^{\ast}(\bm{\beta},\mathbf{B}_{n})=R^{\ast}(\bm{\beta}^{\prime},\mathbf{B}_{n}^{\prime})=$
$\displaystyle\frac{1}{NL}\left[\left\|\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}\left\\{\mathbf{D}_{n}^{-1/2}\left(\mathbf{R}_{n}\circ\mathbf{P}_{n}^{\dagger}\circ\mathbf{Y}_{n}\right)\right\\}-\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}\bm{\beta}^{\prime}\right\|_{F}^{2}+\right.$
$\displaystyle\left.\left\|\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}^{\perp}\left\\{\mathbf{D}_{n}^{-1/2}\left(\mathbf{R}_{n}\circ\mathbf{P}_{n}^{\dagger}\circ\mathbf{Y}_{n}\right)\right\\}-\mathbf{B}_{n}^{\prime}\right\|_{F}^{2}\right].$

It can be seen that $\bm{\beta}^{\ast\prime}$ and
$\mathbf{B}_{n}^{\ast\prime}$ are the unique minimizers of
$R^{\ast}(\bm{\beta}^{\prime},\mathbf{B}_{n}^{\prime})$. Although
$\bm{\beta}^{\ast}$ and $\mathbf{B}_{n}^{\ast}$ cannot be uniquely determined,
we ensure that

$\mathbf{X}_{n}\bm{\beta}^{\ast\prime}+\mathbf{D}_{n}^{1/2}\mathbf{B}_{n}^{\ast\prime}=\mathbf{X}_{n}\bm{\beta}^{\ast}+\mathbf{B}_{n}^{\ast},$

which is sufficient to identify the parameters of interest $\theta_{j}$ for
$j=1,\ldots,L$. Therefore, in what follows, we will focus on estimating
$\bm{\beta}^{\ast\prime}$ and $\mathbf{B}_{n}^{\ast\prime}$.

###  2.3 Estimation of $\bm{\beta}^{\ast\prime}$ and
$\mathbf{B}_{n}^{\ast\prime}$

Because $\mathbf{P}_{n}$ is unknown, we consider a maximum likelihood
estimator $\widehat{\mathbf{P}}_{n}$ of $\mathbf{P}_{n}$ and

$\displaystyle\widehat{R}^{\ast}(\bm{\beta}^{\prime},\mathbf{B}_{n}^{\prime})=$
$\displaystyle\frac{1}{NL}\left[\left\|\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}\left\\{\mathbf{D}_{n}^{-1/2}\left(\mathbf{R}_{n}\circ\widehat{\mathbf{P}}_{n}^{\dagger}\circ\mathbf{Y}_{n}\right)\right\\}-\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}\bm{\beta}^{\prime}\right\|_{F}^{2}+\right.$
$\displaystyle\left.\left\|\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}^{\perp}\left\\{\mathbf{D}_{n}^{-1/2}\left(\mathbf{R}_{n}\circ\widehat{\mathbf{P}}_{n}^{\dagger}\circ\mathbf{Y}_{n}\right)\right\\}-\mathbf{B}_{n}^{\prime}\right\|_{F}^{2}\right],$

where $\widehat{\mathbf{P}}_{n}^{\dagger}$ is the matrix of the estimated
response probabilities. Since $\bm{\beta}^{\prime}$ and
$\mathbf{B}_{n}^{\prime}$ are high-dimensional parameters, a direct
minimization of $\widehat{R}^{\ast}(\bm{\beta},\mathbf{B}_{n})$ would often
result in over-fitting. To avoid such an issue, we incorporate penalty terms
for those two parameters. Specifically, we propose the penalized estimators of
$(\bm{\beta}^{\ast\prime},\mathbf{B}_{n}^{\ast\prime})$ as

$(\widehat{\bm{\beta}}^{\prime},\widehat{\mathbf{B}}_{n}^{\prime})=\operatorname*{arg\,min}_{\bm{\beta}^{\prime},\mathbf{B}_{n}^{\prime}\in\mathcal{N}(\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n})}\widehat{R}^{\ast}(\bm{\beta}^{\prime},\mathbf{B}_{n}^{\prime})+\tau_{1}\left\|\bm{\beta}^{\prime}\right\|_{F}^{2}+\tau_{2}\left\\{\alpha\left\|\mathbf{B}_{n}^{\prime}\right\|_{\ast}+(1-\alpha)\left\|\mathbf{B}_{n}^{\prime}\right\|_{F}^{2}\right\\},$
(8)

where $\|\mathbf{M}\|_{\ast}=\mbox{trace}(\sqrt{\mathbf{M}^{\rm
T}\mathbf{M}})$ is the nuclear norm of a real-valued matrix $\mathbf{M}$, and
$\tau_{1}$, $\tau_{2}>0$ along with $0\leq\alpha\leq 1$ are regularization
parameters. Since $\mathbf{B}_{N}^{\ast}$ is assumed to be low-rank,
$\mathbf{B}_{n}^{\ast}$ is also low-rank and
$\text{rank}(\mathbf{B}_{n}^{\ast\prime})=\text{rank}(\mathbf{B}_{n}^{\ast})$.
Similar to the rank sum norm, the nuclear norm also encourages a low-rank
solution. In matrix completion literature, one reason why people consider the
nuclear norm instead of the rank norm penalty directly is that the
minimization problem with rank norm penalty is NP-hard (Candès and Recht,
2009). The two additional Frobenius norm terms of $\bm{\beta}^{\prime}$ and
$\mathbf{B}_{n}^{\prime}$ are applied to improve finite sample performance
(Zou and Hastie, 2005; Sun and Zhang, 2012; Mao, Chen and Wong, 2019).

To obtain $\widehat{\bm{\beta}}^{\prime}$, it is essentially a solution of a
ridge regression problem, and we have

$\widehat{\bm{\beta}}^{\prime}=\left(\mathbf{X}_{n}^{\rm
T}\mathbf{D}_{n}^{-1}\mathbf{X}_{n}+NL\tau_{1}\mathbf{I}\right)^{-1}\mathbf{X}_{n}^{\rm
T}\mathbf{D}_{n}^{-1}\left(\mathbf{R}_{n}\circ\widehat{\mathbf{P}}_{n}^{\dagger}\circ\mathbf{Y}_{n}\right).$

To obtain $\widehat{\mathbf{B}}_{n}^{\prime}$, following the same argument in
Proposition 2 of Mao, Chen and Wong (2019), we can extend the searching domain
for
$\mathbf{B}_{n}^{\prime}\in\mathcal{N}(\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n})$
in the minimization problem (8) to be
$\mathbf{B}_{n}^{\prime}\in\mathbb{R}^{n\times L}$. This allows us to express
the solution $\widehat{\mathbf{B}}_{n}^{\prime}$ in a closed form. Let
$\mathbf{U}\bm{\Sigma}\mathbf{V}^{\rm T}$ be the SVD of a matrix $\mathbf{M}$,
where $\bm{\Sigma}=\mathrm{diag}(\\{\sigma_{i}\\})$. Define the corresponding
singular value soft-thresholding operator $\mathcal{T}_{c}$ by

$\mathcal{T}_{c}\left(\mathbf{M}\right)=\mathbf{U}\mathrm{diag}(\\{\left(\sigma_{i}-c\right)_{+}\\})\mathbf{V}^{\intercal}$

for any $c\geq 0$, where $x_{+}=\max(x,0)$. It can be shown that the solution
$\widehat{\mathbf{B}}_{n}^{\prime}$ in (8) possesses the following closed
form:

$\widehat{\mathbf{B}}_{n}^{\prime}=\frac{1}{1+\left(1-\alpha\right)NL\tau_{2}}\mathcal{T}_{\alpha
NL\tau_{2}/2}\left[\mathcal{P}_{\mathbf{D}_{n}^{-1/2}\mathbf{X}_{n}}^{\perp}\left\\{\mathbf{D}_{n}^{-1/2}\left(\mathbf{R}_{n}\circ\widehat{\mathbf{P}}_{n}^{\dagger}\circ\mathbf{Y}_{n}\right)\right\\}\right].$

Following the common practice in matrix completion works (Mazumder et al.,
2010; Xu, Jin and Zhou, 2013; Chiang, Hsieh and Dhillon, 2015; Mao, Wong and
Chen, 2019), we obtain tuning parameters $\tau_{1}$, $\tau_{2}$ and $\alpha$
by a $5$-fold cross validation procedure. After obtaining
$(\widehat{\bm{\beta}}^{\prime},\widehat{\mathbf{B}}_{n}^{\prime})$, an
estimator of $\mathbf{A}_{n}$ is given by

$\widehat{\mathbf{A}}_{n}=\mathbf{X}_{n}\widehat{\bm{\beta}}^{\prime}+\mathbf{D}_{n}^{1/2}\widehat{\mathbf{B}}_{n}^{\prime}.$
(9)

###  2.4 Comparison with Hot Deck Imputation and Multiple Imputation

It is worth comparing the proposed matrix completion method with existing
approaches for imputation. Hot deck imputation uses an observed datum as a
“donor” to impute each missing item based on a specific distance using some
fully observed auxiliary information. For hot deck imputation, an underlying
regression model, $f_{j}(\bm{x}_{i})$, is assumed for the item $y_{ij}$.
Therefore, only $\bm{x}_{i}$ is used for imputing $y_{ij}$ but not $y_{ik}$
with $k\neq j$. The multiple imputation (Rubin, 1978) assumes a joint model of
$(\bm{x}_{i},\bm{y}_{i})$ and uses all available variables for imputation.
However, fully parametric modeling is sensitive to model misspecification.

In our approach, the low-rank structure of $\mathbf{A}_{N}$ suggests a general
decomposition of $\mathbf{A}_{N}$ to be
$\mathbf{A}_{N}=\mathbf{U}_{N}\mathbf{V}_{N}^{\rm T}$, where
$\mathbf{U}_{N}\in\mathbb{R}^{N\times r_{\mathbf{A}_{N}}}$ and
$\mathbf{V}_{N}\in\mathbb{R}^{L\times r_{\mathbf{A}_{N}}}$ are two hidden
matrices. Due to the low-rank assumption, we have $r_{\mathbf{A}_{N}}\ll N$
and $r_{\mathbf{A}_{N}}\ll L$. In our column-space-decomposition model, we
enforce part of the hidden matrix $\mathbf{U}_{N}$ to be a fully observed
matrix $\mathbf{X}_{N}\in\mathbb{R}^{N\times d}$ and denote the corresponding
part in $\mathbf{V}_{N}$ to be $\bm{\beta}^{\ast}$, where $\bm{\beta}^{\ast}$
is just a different notation and still totally unknown. Thus, the
decomposition could be written as
$\mathbf{A}_{N}=(\mathbf{X}_{N},\mathbf{U}_{N}^{\ast})(\bm{\beta}^{\ast},\mathbf{V}_{N}^{\ast})^{\rm
T}$ with
$\mathbf{B}_{N}^{\ast}=\mathbf{U}_{N}^{\ast}{\mathbf{V}_{N}^{\ast}}^{\rm T}$.
In a general setting, the only restriction for $\mathbf{U}_{N}^{\ast}$ is
$\text{rank}(\mathbf{X}_{N},\mathbf{U}_{N}^{\ast})=r_{\mathbf{A}_{N}}$, which
means that each column of $\mathbf{U}_{N}^{\ast}$ cannot be fully expressed by
the columns in $\mathbf{X}_{N}$. However, it still allows for
$\text{cor}(\mathbf{X}_{N},\mathbf{U}_{N}^{\ast})\neq\mathbf{0}$. Then, it is
difficulty to identify the hidden matrix $\mathbf{U}_{N}^{\ast}$ under the
general setting. Thus, we restrict the column space of $\mathbf{U}_{N}^{\ast}$
to be orthogonal to the column space of $\mathbf{X}_{N}$. Fortunately, the
number of covariates $d$ is usually fixed and $d\ll r_{\mathbf{A}_{N}}$, which
means that we would not lose too much freedom for $\mathbf{U}_{N}^{\ast}$.

###  2.5 Estimation of $\theta_{j}$

After imputation, it may be natural to estimate $\theta_{j}$ by the Horvitz-
Thompson estimator (5) applied to the imputed dataset. However, it is well
known that the estimated low-rank matrix $\widehat{\mathbf{B}}_{n}^{\prime}$
is biased when $n$ is finite (Mazumder et al., 2010; Foucart et al., 2017;
Carpentier and Kim, 2018; Chen et al., 2019). Therefore, the resulting
imputation estimator is biased. Researchers have proposed different procedures
to alleviate or eliminate the bias. Mazumder et al. (2010) suggested a post-
processing step by re-estimating the estimated singular values without any
theoretical guarantee. Foucart et al. (2017) proposed an algorithm based on
projection onto the max-norm ball to de-bias the estimator under non-uniform
and deterministic sampling patterns. Carpentier and Kim (2018) considered an
estimator using an iterative hard thresholding method and showed that the
entry-wise bias is small when the sampling design is Gaussian. More recently,
Chen et al. (2019) developed a de-biasing procedure using a similar idea to
de-biasing LASSO estimators and showed nearly optimal properties for the
resulting estimator. Despite these advances in literature, the scenarios
considered above are restricted to deterministic sampling, Gaussian sampling
or missing completely at random (MCAR), which are not applicable in our
setting.

We use a simple strategy borrowing the idea from the doubly robust estimation
literature (Robins, Rotnitzky and Zhao, 1994; Bang and Robins, 2005; Cao,
Tsiatis and Davidian, 2009) and consider a doubly robust estimator of
$\theta_{j}$ as

$\widehat{\theta}_{j,DR}=\frac{1}{N}\sum_{i=1}^{N}\frac{I_{i}}{\pi_{i}}\left\\{\frac{r_{ij}(y_{ij}-\widehat{a}_{ij})}{\widehat{p}_{ij}}+\widehat{a}_{ij}\right\\},$
(10)

where $\widehat{p}_{ij}$ and $\widehat{a}_{ij}$ are the $(i,j)$th element of
$\widehat{P}_{n}$ and $\widehat{A}_{n}$, respectively. It can be shown that

$\widehat{\theta}_{j,DR}=\frac{1}{N}\sum_{i=1}^{N}\frac{I_{i}}{\pi_{i}}\left\\{\frac{r_{ij}(y_{ij}-\widehat{a}_{ij})}{{p}_{ij}}+\widehat{a}_{ij}\right\\}+o_{P}(1)$

when the response model (3) is correctly specified, so
$\widehat{\theta}_{j,DR}$ is asymptotically unbiased for $\theta_{j}$ for this
case.

##  3 Asymptotic Properties

In this section, we first study the asymptotic properties of the estimator
$\widehat{\mathbf{A}}_{n}$ in (9) under the logistic regression model (3).
Further, we establish the average convergence rate of
$\widehat{\theta}_{j,DR}-\theta_{j}$ for $j=1,\ldots,L$.

For asymptotic inference, we follow the framework of Isaki and Fuller (1982)
and assume that both the population size $N$ and the sample size $n$ go to
infinity. Let $\|\mathbf{M}\|=\sigma_{\max}(\mathbf{M})$ and
$\|\mathbf{M}\|_{\infty}=\max_{i,j}|m_{ij}|$ be the spectral and the maximum
norms of a matrix $\mathbf{M}$, respectively. We use the symbol “$\asymp$” to
represent the asymptotic equivalence in order, that is, $a_{n}\asymp b_{n}$ is
equivalent to $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$.

The technical conditions needed for our analysis are given as follows.

  1. C1

(a) The random errors $\\{\epsilon_{ij}\\}_{i,j=1}^{N,L}$ in (2) are
independently distributed random variables such that $E(\epsilon_{ij})=0$ and
$E(\epsilon^{2}_{ij})=\sigma_{ij}^{2}<\infty$ for all $i,j$. (b) For some
finite positive constants $c_{\sigma}$ and $\eta$,
$\underset{i,j}{\max}\,E|\epsilon_{ij}|^{l}\leq\frac{1}{2}l!c_{\sigma}^{2}\eta^{l-2}$
for any positive integer $l\geq 2$.

  2. C2

The inclusion probability satisfies $\pi_{i}\asymp nN^{-1}$ for
$i=1,\ldots,N$.

  3. C3

The population design matrix $\mathbf{X}_{N}$ is of size $N\times d$ such that
$N>d$. Moreover, there exists a positive constant $a_{x}$ such that
$\|\mathbf{X}_{N}\|_{\infty}\leq a_{x}$ and $\mathbf{X}_{N}^{\rm
T}\mathbf{D}_{N}\mathbf{X}_{N}$ is invertible, where $\mathbf{D}_{N}$ is a
diagonal matrix with $\pi_{i}$ as its $(i,i)$th entry. Furthermore, there
exists a symmetric matrix $\mathbf{S}_{\mathbf{X}}$ with
$\sigma_{\min}(\mathbf{S}_{\mathbf{X}})\asymp
1\asymp\|\mathbf{S}_{\mathbf{X}}\|$ such that $n_{0}^{-1}\mathbf{X}_{N}^{\rm
T}\mathbf{D}_{N}\mathbf{X}_{N}\to\mathbf{S}_{\mathbf{X}}$ as $N\to\infty$,
where $n_{0}=\sum_{i=1}^{N}\pi_{i}$ is the expected sample.

  4. C4

There exists a positive constant $a$ such that
$\max\\{\|\mathbf{X}_{N}\bm{\beta}^{\ast}\|_{\infty},\|\mathbf{A}_{N}\|_{\infty}\\}\leq
a$.

  5. C5

The indicators of observed entries $\\{r_{ij}\\}_{i,j=1}^{N,L}$ are mutually
independent, $r_{ij}\sim\text{Bern}(p_{ij})$ for $p_{ij}\in(0,1)$ and are
independent of $\\{\epsilon_{ij}\\}_{i,j=1}^{N,L}$ given $\mathbf{X}_{N}$.
Furthermore, for $i=1,\dots,N$ and $j=1,\dots,L$,
$\Pr(r_{ij}=1|\bm{x}_{i},y_{ij})=\Pr(r_{ij}=1|\bm{x}_{i})$ follows the
logistic regression model (3).

  6. C6

There exists a lower bound $p_{\min}\in(0,1)$ such that
$\underset{i,j}{\min}\\{p_{ij}\\}\geq p_{\min}>0$, where $p_{\min}$ is allowed
to depend on $n$ and $L$. The number of questions $L\leq n$.

Condition C1(a) is a common regularity condition for the measurement errors in
$\bm{\epsilon}_{N}$, and C1(b) is the Bernstein condition (Koltchinskii et
al., 2011). Condition C2 is widely used in survey sampling and regulates the
inclusion probabilities of a sampling design (Fuller, 2009). To illustrate
ideas, we consider Poisson sampling in this section, and our discussion
applies to other sampling designs such as simple random sampling and
probability-proportional-to-size sampling. In Condition C3, the requirement
$N>d$ is easily met as the number of questions in a survey is usually fixed,
and the population size is often larger than the number of questions. As the
dimension of $n_{0}^{-1}\mathbf{X}_{N}^{\rm T}\mathbf{D}_{N}\mathbf{X}_{N}$ is
fixed at $d\times d$, it is mild to assume $\mathbf{X}_{N}^{\rm
T}\mathbf{D}_{N}\mathbf{X}_{N}$ to be invertible, and there exists a symmetric
matrix $\mathbf{S}_{\mathbf{X}}$ as the probability limit of
$n_{0}^{-1}\mathbf{X}_{N}^{\rm T}\mathbf{D}_{N}\mathbf{X}_{N}$. Furthermore,
the sample size is often larger than the number of questions, that is, $n>d$,
and it is not hard to show that together with Condition C2, the probability
limit of $n^{-1}\mathbf{X}_{n}^{\rm T}\mathbf{X}_{n}$ is also
$\mathbf{S}_{\mathbf{X}}$ under Poisson sampling; see the Supplementary
Materials (Mao, Wang and Yang, 20xx) for details. The order of
$\sigma_{\min}(\mathbf{S}_{\mathbf{X}})$ and $\|\mathbf{S}_{\mathbf{X}}\|$
equals to $1$ is due to $\|\mathbf{X}_{N}\|_{\infty}<\infty$. Condition C4 is
also standard in the matrix completion literature (Koltchinskii et al., 2011;
Negahban and Wainwright, 2012; Cai and Zhou, 2016). Especially, it is
reasonable to assume all the responses are bounded in survey sampling.
Condition C5 describes the independent Bernoulli model for the response
indicator of observing $y_{ij}$, where the probability of observation $p_{ij}$
follows the logistic model (3). In Condition C6, the lower bound $p_{\min}$ is
allowed to go to $0$ with $n$ and $L$ growing. This condition is more general
than we need for a typical survey, and $p_{\min}\asymp 1$ suffices. Typically,
the number of questions $L$ grows slower than the number of participants $n$
in survey sampling. Thus, the assumption that $L\leq n$ is quite mild.

For any $\delta_{\sigma}>0$, some positive constants $C_{d}$, $C_{g}$, $C$ and
$t\in(d+3,+\infty)$, define

$\Delta\left(\delta_{\sigma},t\right)=\max\left\\{N^{1/2}n^{-1}L^{-1}\log^{1/2}\left(n\right)p_{\min}^{-1/2},N^{1/2}n^{-5/4}L^{-1/4}\log^{1/2}\left(L\right)\log^{\delta_{\sigma}/4}\left(n\right)t^{1/2}p_{\min}^{-3/2}\right\\},$
(11)

and
$\eta_{n,L}(\delta_{\sigma},t)=4/(n+L)+4C_{d}t\exp\\{-t/2\\}+4/L+C\log^{-\delta_{\sigma}}(n)$.
We can verify that
$\lim_{t\to\infty}\\{\lim_{n,L\to\infty}\eta_{n,L}(\delta_{\sigma},t)\\}=0$.
Once we have $n^{1/2}L^{-3/2}\log\left(n\right)p_{\min}^{2}\geq(d+3)$, by
choosing $t$ such that

$d+3<t<n^{1/2}L^{-3/2}\log\left(n\right)p_{\min}^{2},$ (12)

we can show $\underset{t}{\sup}\Delta(\delta_{\sigma},t)\asymp
N^{1/2}n^{-1}L^{-1}\log^{1/2}(n)p_{\min}^{-1/2}$, which is denoted by
$\Delta(\delta_{\sigma})$. Here, the requirement
$n^{1/2}L^{-3/2}\log\left(n\right)p_{\min}^{2}\geq(d+3)$ is easy to fulfill
once $n$ large enough.

###### Theorem 1

Assume Conditions C1-C6 and Poisson sampling,
$p_{\min}^{-1}=O(L\log^{-1}(n+L))$ and the logistic model (3) hold. Choose $t$
as (12), $\tau_{1}\asymp
N^{-1}nL^{-1}{\log^{-1/2}(n)}\Delta(\delta_{\sigma})$,
$1-\alpha\asymp(nL)^{-1}$, $\tau_{2}\asymp
p_{\min}^{-3/2}N^{-1}n^{1/4}L^{-1/4}\log^{1/2}(L)\log^{\delta_{\sigma}/3}(n)$
in (8) for any $\delta_{\sigma}>0$. Then, for some positive constant $C_{1}$
and $C_{2}$, with probability at least $1-\eta_{n,L}(\delta_{\sigma},t)$, we
have

$\displaystyle\quad\frac{1}{mL}\left\|\widehat{\bm{\beta}}^{\prime}-\bm{\beta}^{\ast\prime}\right\|_{F}^{2}\leq
C_{1}r_{\mathbf{B}_{N}}L^{-1}\log\left(n\right)p_{\min}^{-1}\quad\text{and}$
$\displaystyle\quad\frac{1}{nL}\left\|\widehat{\mathbf{B}}_{n}^{\prime}-\mathbf{B}_{n}^{\ast\prime}\right\|_{F}^{2}\leq
C_{2}r_{\mathbf{B}_{N}}Nn^{-1}L^{-1}\log\left(n\right)p_{\min}^{-1}.$

A proof of Theorem 1 is given in the Supplementary Materials (Mao et al.,
20xx). As
$\lim_{t\to\infty}\\{\lim_{n,L\to\infty}\eta_{n,L}(\delta_{\sigma},t)\\}=0$,
Theorem 1 implies that
$(mL)^{-1}\|\widehat{\bm{\beta}}^{\prime}-\bm{\beta}^{\ast\prime}\|_{F}^{2}=O_{p}\\{r_{\mathbf{B}_{N}}L^{-1}\log(n)p_{\min}^{-1}\\}$
and
$(nL)^{-1}\|\widehat{\mathbf{B}}_{n}^{\prime}-\mathbf{B}_{n}^{\ast\prime}\|_{F}^{2}=O_{p}\\{r_{\mathbf{B}_{N}}Nn^{-1}L^{-1}\log(n)p_{\min}^{-1}\\}$.

As we pointed out in Section 2.2, even with the knowledge of
$(\bm{\beta}^{\ast\prime},\mathbf{B}_{n}^{\ast\prime})$, we cannot recover
$(\bm{\beta}^{\ast},\mathbf{B}_{n}^{\ast})$ exactly. Fortunately, we have

$\widehat{\mathbf{A}}_{n}=\mathbf{X}_{n}\widehat{\bm{\beta}}^{\prime}+\mathbf{D}_{n}^{1/2}\widehat{\mathbf{B}}_{n}^{\prime},$

which enables us to derive the asymptotic bound for
$(nL)^{-1}\|\widehat{\mathbf{A}}_{n}-\mathbf{A}_{n}\|_{F}^{2}$ given in the
following theorem.

###### Theorem 2

Assume that the same conditions in Theorem 1 hold. For a positive constant
$C_{3}$, with probability at least $1-\eta_{n,L}(\delta_{\sigma},t)$, we have

$\displaystyle\frac{1}{nL}\left\|\widehat{\mathbf{A}}_{n}-\mathbf{A}_{n}\right\|_{F}^{2}\leq
C_{3}r_{\mathbf{B}_{N}}L^{-1}\log\left(n\right)p_{\min}^{-1}.$

A brief proof of Theorem 2 can be found in the Supplementary Materials. The
term $(nL)^{-1}\|\widehat{\mathbf{A}}_{n}-\mathbf{A}_{n}\|_{F}^{2}$ has the
same order with upper bound of
$(mL)^{-1}\|\widehat{\bm{\beta}}^{\prime}-\bm{\beta}^{\ast\prime}\|_{F}^{2}$.
To ensure the convergence of
$(nL)^{-1}\|\widehat{\mathbf{A}}_{n}-\mathbf{A}_{n}\|_{F}^{2}$, we only
require that $n=O\\{\exp(r_{\mathbf{B}_{N}}^{-1}Lp_{\min})\\}$ which is quite
mild. In survey sampling, it is reasonable to assume that $p_{\min}\asymp 1$,
especially when the participants are awarded. Thus, the assumption that
$p_{\min}^{-1}=O(L\log^{-1}(n+L))$ is easy to fulfill once $L$ large enough.
It can be shown that the convergence rate for
$(nL)^{-1}\|\widehat{\mathbf{A}}_{n}-\mathbf{A}_{n}\|_{F}^{2}$ can be
simplified to $r_{\mathbf{B}_{N}}L^{-1}\log(n)$ if $p_{\min}\asymp 1$. As we
have discussed in Section 2.4, the proposed method achieves robustness against
model misspecification.

The following theorem provides the average convergence rate of
$\widehat{\theta}_{j,DR}$ for $j=1,\ldots,L$.

###### Theorem 3

Assume that the same conditions in Theorem 1 hold and $p_{\min}\asymp 1$.
Then, we have

$\frac{1}{L}\sum_{j=1}^{L}(\widehat{\theta}_{j,DR}-\theta_{j})^{2}=O_{p}\\{r_{\mathbf{B}_{N}}L^{-1}\log\left(n\right)\\}.$

A proof for Theorem 3 is given in the Supplementary Materials. By Theorem 3,
the mean squared difference between $\widehat{\theta}_{j,DR}$ and $\theta_{j}$
among the $L$ questions is bounded by
$O_{p}\\{r_{\mathbf{B}_{N}}L^{-1}\log\left(n\right)\\}$. To ensure the
convergence of $L^{-1}\sum_{j=1}^{L}(\widehat{\theta}_{j,DR}-\theta_{j})^{2}$,
similarly with before, we only require that
$n=O\\{\exp(r_{\mathbf{B}_{N}}^{-1}L)\\}$ which is quite mild.

##  4 Simulation

We use (1) and (2) to generate a finite population $U_{N}$, where elements in
$\mathbf{X}_{N}$ and $\bm{\beta}^{\ast}$ are generated by
$\mathcal{N}(0.5,1^{2})$,
$\mathbf{B}_{N}^{\ast}=\mathcal{P}_{\mathbf{X}_{N}}^{\perp}\mathbf{B}_{L}\mathbf{B}_{R}$,
$\mathbf{B}_{L}$ is an $N\times k$ matrix, $\mathbf{B}_{R}$ is a $k\times L$
matrix, elements of $\mathbf{B}_{L}$ and $\mathbf{B}_{R}$ are generated by
$\mathcal{N}(1,3^{2})$, elements of $\bm{\epsilon}_{N}$ are generated such
that the signal-noise ratio is 2, $N=10\,000$ is the population size, $L=500$
is the number of questions in the survey, $d=20$ is the rank of
$\mathbf{X}_{N}$ and $\bm{\beta}^{\ast}$, and $k=10$ is the rank of
$\mathbf{B}_{N}^{\ast}$. From the generated finite population $U_{N}$, the
following sampling designs are considered:

  1. I

Poisson sampling with inclusion probability
$\pi_{i}=ns_{i}(\sum_{i=1}^{N}s_{i})^{-1}$, where $s_{i}>0$ is a size measure
of the $i$th subject, and the generation of $s_{i}$ is discussed later.
Specifically, for $i=1,\ldots,N$, a sampling indicator $I_{i}$ is generated by
a Bernoulli distribution with success probability $\pi_{i}$.

  2. II

Simple random sampling with sample size $n$.

  3. III

Probability-proportional-to-size sampling with size measure $s_{i}$. That is,
a sample of size $n$ is selected independently from the finite population
$U_{N}$ with replacement, and the selection probability of the $i$th subject
is proportional to its size measure $s_{i}$.

We consider two scenarios for the sampling procedure. One is informative
sampling with $s_{i}=7^{-1}\sum_{j=1}^{7}y_{ij}-m_{s}+1$, where
$m_{s}=\min\\{7^{-1}\sum_{j=1}^{7}y_{ij}:i=1,\ldots,N\\}$. The other is
noninformative sampling with $s_{i}=d^{-1}\sum_{j=1}^{d}x_{ij}+e_{i}+1$, where
$e_{i}\sim\mbox{Ex}(1)$, and $\mbox{Ex}(\lambda)$ is an exponential
distribution with rate parameter $\lambda$. Two different sample sizes are
considered, $n=200$ and $n=500$, and the following estimation methods are
compared:

  1. I

Hot deck imputation. For each item with $r_{ij}=0$, we use $y_{kj}$ as the
imputed value, where $\bm{x}_{k}$ is nearest to $\bm{x}_{j}$ among
$\\{\bm{x}_{l}:r_{lj}=1\\}$. Treating the imputed values as observed ones, we
estimate $\theta_{j}$ by (5).

  2. II

Multiple imputation. We adopt the multivariate imputation by chained equations
(MICE) by van Buuren and Groothuis-Oudshoorn (2011). MICE fully specifies the
conditional distribution for the missing data and uses a posterior predictive
distribution to generate imputed values for the nonresponse items; check van
Buuren and Groothuis-Oudshoorn (2011) for details. However, it is impossible
for MICE to impute all missing responses in $\mathbf{Y}_{n}$ at the same time
due to the computational issues. For comparison, we only use the first 20
items of $\mathbf{Y}_{n}$ to specify the conditional distribution for MICE and
generate imputed values for the corresponding nonresponses. Then, we can use
(5) to estimate $\theta_{j}$.

  3. III

Inverse probability method. For $j=1,\ldots,L$, a logistic regression model
(3) is fitted. Then, $\theta_{j}$ is estimated by

$\widehat{\theta}_{j,IPM}=N^{-1}\sum_{i=1}^{n}r_{ij}\widehat{p}_{ij}^{-1}y_{ij}.$

  4. IV

Doubly robust estimator using linear regression model. For $j=1,\ldots,L$,
consider the following linear regression model:

$y_{ij}=\phi_{0j}+\bm{x}_{i}^{\rm T}\bm{\phi}_{1j},$ (13)

and the parameters in (13) are estimated by

$(\widehat{\phi}_{0j},\widehat{\bm{\phi}}_{1j})=\operatorname*{arg\,min}_{(\phi_{0j},\bm{\phi}_{1j})}\sum_{i=1}^{n}\frac{r_{ij}}{\pi_{i}\widehat{p}_{ij}}(y_{ij}-\phi_{0j}-\bm{x}_{i}^{\rm
T}\bm{\phi}_{1j})^{2}.$

Then, we can use a doubly robust estimator based on the linear model (13) to
estimate $\theta_{j}$.

  5. V

Doubly robust estimator using naive imputation. We use the naive imputation
method (Mazumder et al., 2010) by assuming MCAR to generate the imputed
values, and use the doubly robust estimator to estimate $\theta_{j}$.

  6. VI

Doubly robust estimator using the proposed method in (10).

For comparison, we also consider the Horvitz-Thompson estimator in (5) using
the fully observed data.

We conduct 1 000 Monte Carlo simulations. Table 1 shows the Monte Carlo bias
and standard error for the first five items under informative probability-
proportional-to-size sampling with sample size $n=500$. Specifically, the
Monte Carlo bias and standard error for the $j$th question are obtained by

$\mbox{Bias}_{j}=\widehat{\theta}_{j}^{(m)}-\theta_{j}\qquad\text{and}\qquad\mbox{SE}_{j}=\frac{1}{1\,000}\sum_{m=1}^{1\,000}(\widehat{\theta}_{j}^{(m)}-\widehat{\theta}_{j})^{2},$

respectively, where
$\widehat{\theta}_{j}={1\,000^{-1}}\sum_{m=1}^{1\,000}\widehat{\theta}_{j}^{(m)}$,
and $\widehat{\theta}_{j}^{(m)}$ is an estimator from a specific estimation
method in the $m$th Monte Carlo simulation. The standard error for the hot
deck imputation is much larger compared with other methods. The bias of the
multiple imputation is larger than the inverse probability method and doubly
robust estimators since the model is misspecified. Besides, the multiple
imputation method is not preferable due to the computation complexity,
especially when the number of items is large. Two doubly robust estimators
using naive imputation and the proposed method have smaller variability than
others, and the bias for the doubly robust estimator using the proposed method
is smaller. Compared with the Horvitz-Thompson estimator using the fully
observed data, the variance of the doubly robust estimator using the proposed
method is larger.

Table 1: Monte Carlo bias (Bias) and standard error (SE) for the first five items under informative probability-proportional-to-size sampling with sample size $n=500$. “HDI” is the hot deck imputation, “IPM” is the inverse probability method, “DRLR” is the doubly robust estimator using linear regression, “DRNI” is the doubly robust estimator using naive imputation, “DRMC” is the doubly robust estimator using the proposed method, and “Full” is the Horvitz-Thompson estimator using fully observed data. Method | Stat. | Items
---|---|---
I | II | III | IV | V
HDI | Bias | -0.13 | 1.40 | 0.36 | 1.35 | 1.16
SE | 8.03 | 12.20 | 13.19 | 12.75 | 9.54
MI | Bias | 0.29 | 0.64 | 0.30 | 0.69 | 0.56
SE | 1.03 | 1.53 | 1.60 | 1.62 | 1.18
IPM | Bias | -0.03 | -0.02 | -0.01 | 0.05 | 0.14
SE | 1.07 | 1.74 | 1.81 | 1.86 | 1.26
DRLR | Bias | 0.00 | 0.11 | 0.01 | 0.16 | 0.21
SE | 1.07 | 1.76 | 1.82 | 1.88 | 1.26
DRNI | Bias | -0.26 | -0.25 | -0.53 | -0.51 | -0.38
SE | 0.94 | 1.43 | 1.47 | 1.50 | 1.06
DRMC | Bias | -0.16 | 0.02 | -0.21 | -0.11 | -0.04
SE | 0.94 | 1.38 | 1.43 | 1.48 | 1.03
Full | Bias | -0.04 | -0.01 | -0.00 | 0.01 | 0.06
SE | 0.78 | 1.27 | 1.33 | 1.36 | 0.92

Next, we compare different estimation methods by the Monte Carlo mean squared
error (MSE)

$\mbox{MSE}_{j}=\frac{1}{1\,000}\sum_{m=1}^{1\,000}(\widehat{\theta}_{j}^{(m)}-\theta_{j})^{2}\quad(j=1,\ldots,L).$

The result for multiple imputation is omitted due to the computational issue.
Table 2 summarizes the mean and standard error of MSEs for different
questions. From Table 2, we can conclude that the mean MSE and its standard
error of the doubly robust estimator using the proposed method are smallest
among alternatives for all scenarios. Besides, the average MSE and its
standard error of doubly robust estimator using the proposed method are
slightly larger than the Horvitz-Thompson estimator using fully observed data.

Table 2: Summary of MSE for different estimation methods. “NIF” shows the results under noninformative sampling, and “IF” shows those under informative sampling. “POI” for the Poisson sampling, “SRS” for the simple random sampling, and “PPS” stands for the probability-proportional-to-size sampling. “Mean” and “SE” are the mean and standard error of the MSEs for $L=500$ items, “HDI” is the hot deck imputation, “IPM” is the inverse probability method, “DRLR” is the doubly robust estimator using linear regression, “DRNI” is the doubly robust estimator using naive imputation, “DRMC” is the doubly robust estimator using the proposed matrix completion method, and “Full” is the Horvitz-Thompson estimator using fully observed data.  | Design | Sample size | Stat. | Estimation methods
---|---|---|---|---
| HDI | IPM | DRLR | DRNI | DRMC | Full
NIF | POI | $n=200$ | Mean | 17.48 | 11.70 | 12.48 | 8.36 | 7.65 | 6.26
SE | 8.41 | 4.30 | 4.70 | 2.94 | 2.48 | 2.26
$n=500$ | Mean | 15.34 | 4.54 | 4.63 | 3.45 | 2.92 | 2.45
SE | 7.69 | 1.70 | 1.75 | 1.39 | 0.99 | 0.89
SRS | $n=200$ | Mean | 16.91 | 9.97 | 10.91 | 7.28 | 6.36 | 5.11
SE | 8.26 | 3.81 | 4.00 | 2.72 | 2.05 | 1.91
$n=500$ | Mean | 15.12 | 3.82 | 3.96 | 3.01 | 2.44 | 1.96
SE | 7.79 | 1.46 | 1.49 | 1.23 | 0.81 | 0.76
PPS | $n=200$ | Mean | 17.10 | 11.16 | 12.04 | 7.98 | 7.10 | 5.82
SE | 8.32 | 4.29 | 4.56 | 2.93 | 2.35 | 2.10
$n=500$ | Mean | 15.32 | 4.47 | 4.57 | 3.37 | 2.85 | 2.37
SE | 7.84 | 1.64 | 1.69 | 1.34 | 0.91 | 0.88
IF | POI | $n=200$ | Mean | 17.25 | 11.23 | 12.07 | 8.11 | 6.98 | 5.97
SE | 8.17 | 4.26 | 4.51 | 2.98 | 2.24 | 2.14
$n=500$ | Mean | 15.34 | 4.35 | 4.43 | 3.32 | 2.70 | 2.32
SE | 7.70 | 1.68 | 1.70 | 1.38 | 0.84 | 0.81
SRS | $n=200$ | Mean | 16.79 | 9.86 | 10.80 | 7.32 | 6.40 | 5.11
SE | 8.32 | 3.74 | 4.10 | 2.83 | 2.11 | 1.90
$n=500$ | Mean | 14.99 | 3.83 | 3.91 | 2.96 | 2.41 | 1.97
SE | 7.74 | 1.42 | 1.44 | 1.18 | 0.78 | 0.73
PPS | $n=200$ | Mean | 16.91 | 10.64 | 11.65 | 7.66 | 6.54 | 5.45
SE | 8.01 | 4.07 | 4.26 | 2.64 | 2.08 | 2.02
$n=500$ | Mean | 15.20 | 4.29 | 4.38 | 3.19 | 2.58 | 2.22
SE | 7.70 | 1.60 | 1.61 | 1.24 | 0.82 | 0.81

##  5 Application

The NHANES 2015-2016 Questionnaire Data is used as an application for the
proposed method. Conducted by the National Center for Health Statistics, the
NHANES is a unique survey combining interviews and physical examinations to
study the health and nutritional status of adults and children in the United
States. Data are released in a two-year cycle. The sample size is
approximately 5 000, and the participants are nationally representative. The
sampling design for NHANES aims at reliable estimation for population
subgroups formed by age, sex, income status and origins. Specifically, the
Questionnaire Data contains family-level information including food security
status as well as individual level information including dietary behavior and
alcohol use.

In this section, we are interested in estimating the population mean of
alcohol usage, blood pressure and cholesterol, diet behavior and nutrition,
diabetes status, mental health, income status, and sleep disorders based on
the newly released NHANES 2015-2016 Questionnaire Data. There are about 39
blocks of questions in this dataset. Each block contains several relevant
questions, and the number of questions in our analysis ranges from 4 to 17.
There are $n=5\,735$ eligible subjects involved and 146 items including 45
demographic questions. Among the demographic items, there are 16 fully
observed items including age, gender and race-ethnicity, and they are used as
the covariates $\mathbf{X}_{n}$ in (9). In addition, the sampling weight for
each subject is available. For the questions in our study, the average and
standard error of the missing rates are 0.33 and 0.37, respectively.

For estimating the population mean of each question, we consider those
estimation methods in Section 4, and the covariates are standardized. Since
the population size is unavailable, we use $\widehat{N}=\sum_{i=1}^{n}w_{i}$
instead, where $w_{i}$ is the sampling weight of the $i$th subject
incorporating the sampling design as well as calibration (Fuller, 2009). For
the multiple imputation, we only impute the missing values for the first 20
items due to the computational issue.

Table 3 shows missing rates and estimation results for six randomly selected
items grouped by the missing rate. There are two items with low missing rates
0.08 and 0.09, two with middle missing rates 0.26 and 0.29, and two with high
missing rates 0.65 and 0.68. Besides, there are three items are among the
first 20 used for the multiple imputation. Thus, the selected questions are
representative. When the missing rate is low, estimators are similar for
different methods. As missing rate increases, estimators for the multiple
imputation, hot deck imputation and the double robust estimator using naive
matrix completion are different from those for the inverse probability method
and double robust estimators using linear regression and the proposed method.
When the missing rate is large, say around $.65$, the double robust estimator
using linear regression differs from those for inverse probability method and
the proposed method. Noting that all estimators are unbiased if the response
model is corrected specified; however, the doubly robust estimator with matrix
completion provide the most accurate estimation when all questions are of
interest.

Table 3: Estimation results for six questions. “I” for “Family has savings more than $20,000” “II” stands for “Had at least 12 alcohol drinks/1 yr?”, “III” for “How often drink alcohol over past 12 mos?”, “IV” for “How often drank milk age 5-12?”, “V” for “Told had high blood pressure - 2+ times?”, and “VI” for “Receive community/Government meals delivered?”  Items | Missing rate |  | Estimation methods
---|---|---|---
| MI | HDI | IPM | DRLR | DRNI | DRMC
I | 0.08 |  | - | 1.57 | 1.60 | 1.60 | 1.57 | 1.60
II | 0.09 |  | 1.26 | 1.23 | 1.26 | 1.26 | 1.26 | 1.26
III | 0.26 |  | 3.26 | 2.62 | 3.02 | 3.02 | 2.44 | 3.02
IV | 0.29 |  | - | 2.81 | 2.75 | 2.75 | 2.39 | 2.75
V | 0.65 |  | 1.28 | 1.06 | 1.20 | 1.29 | 0.62 | 1.14
VI | 0.68 |  | - | 1.99 | 0.54 | 3.42 | 0.67 | 0.30

##  6 Concluding Remarks

We have proposed a new imputation method for survey sampling by assuming a
low-rank structure and incorporating fully observed auxiliary information.
Asymptotic properties of the proposed method are investigated. One advantage
of the proposed method is that we can impute the whole survey questionnaire at
the same time. A simulation study demonstrates that the proposed method is
more accurate than some commonly used alternatives, including inverse
probability method and multiple imputation, for estimating all items.

Our framework can also be extended in the following directions. First, we have
considered missingness at random; however, in some situations, the missingness
of $y_{ij}$ may depend on its own value, leading to missingness not at random
(Rubin, 1976); that is, $y_{ij}$ is also involved in the response probability
(3). In this case, we will consider the instrumental variable approach (Wang,
Shao and Kim, 2014; Yang, Wang and Ding, 2019) or stringent parametric model
assumptions (Tang, Little and Raghunathan, 2003; Chang and Kott, 2008; Kim and
Yu, 2011) for identification and estimation. Second, even though we have
proposed an efficient estimator using matrix completion and derived the
asymptotic bounds, its asymptotic distribution is not completely developed,
which will be our future work. Third, because causal inference of treatment
effects can be viewed as a missing data problem, it is intriguing to develop
matrix completion to deal with a partially observed confounder matrix, which
is ubiquitous in practice but has received little attention in the literature
(Yang et al., 2019).

## References

  * (1)
  * Andridge and Little (2010) Andridge, R. R. and Little, R. J. (2010).  A review of hot deck imputation for survey non-response, Int. Stat. Rev. 78(1): 40–64.
  * Bang and Robins (2005) Bang, H. and Robins, J. M. (2005).  Doubly robust estimation in missing data and causal inference models, Biometrics 61: 962–973.
  * Cai and Zhou (2016) Cai, T. T. and Zhou, W.-X. (2016).  Matrix completion via max-norm constrained optimization, Electron. J. Stat. 10(1): 1493–1525.
  * Candès and Recht (2009) Candès, E. J. and Recht, B. (2009).  Exact matrix completion via convex optimization, Found. Comput. Math. 9(6): 717–772.
  * Cao et al. (2009) Cao, W., Tsiatis, A. A. and Davidian, M. (2009).  Improving efficiency and robustness of the doubly robust estimator for a population mean with incomplete data, Biometrika 96: 723–734.
  * Carpentier and Kim (2018) Carpentier, A. and Kim, A. K. (2018).  An iterative hard thresholding estimator for low rank matrix recovery with explicit limiting distribution, Statist. Sinica 28: 1371–1393.
  * Chang and Kott (2008) Chang, T. and Kott, P. S. (2008).  Using calibration weighting to adjust for nonresponse under a plausible model, Biometrika 95: 555–571.
  * Chen and Shao (2000) Chen, J. and Shao, J. (2000).  Nearest neighbor imputation for survey data, J. Off. Stat. 16: 113–131.
  * Chen et al. (2019) Chen, Y., Fan, J., Ma, C. and Yan, Y. (2019).  Inference and uncertainty quantification for noisy matrix completion, arXiv preprint arXiv:1906.04159 .
  * Chiang et al. (2015) Chiang, K.-Y., Hsieh, C.-J. and Dhillon, I. S. (2015).  Matrix completion with noisy side information, Adv. Neural Inf. Process. Syst., Vol. 28, pp. 3447–3455.
  * Clogg et al. (1991) Clogg, C. C., Rubin, D. B., Schenker, N., Schultz, B. and Weidman, L. (1991).  Multiple imputation of industry and occupation codes in census public-use samples using Bayesian logistic regression, J. Amer. Statist. Assoc. 86: 68–78.
  * Davenport and Romberg (2016) Davenport, M. A. and Romberg, J. (2016).  An overview of low-rank matrix recovery from incomplete observations, IEEE J. Sel. Topics Signal Process. 10(4): 608–622.
  * Fay (1991) Fay, R. E. (1991).  A design-based perspective on missing data variance, US Census Bureau [custodian].
  * Fay (1992) Fay, R. E. (1992).  When are inferences from multiple imputation valid?, Proceedings of the Survey Research Methods Section of the American Statistical Association, American Statistical Association, pp. 227–232.
  * Foucart et al. (2017) Foucart, S., Needell, D., Plan, Y. and Wootters, M. (2017).  De-biasing low-rank projection for matrix completion, Wavelets and Sparsity XVII, Vol. 10394, International Society for Optics and Photonics, p. 1039417.
  * Fuller (2009) Fuller, W. A. (2009).  Sampling Statistics, Wiley, Hoboken, NJ.
  * Fuller and Kim (2005) Fuller, W. A. and Kim, J. K. (2005).  Hot deck imputation for the response model, Surv. Methodol. 31: 139.
  * Haziza and Rao (2006) Haziza, D. and Rao, J. N. (2006).  A nonresponse model approach to inference under imputation for missing survey data, Surv. Methodol. 32(1): 53.
  * Horvitz and Thompson (1952) Horvitz, D. G. and Thompson, D. J. (1952).  A generalization of sampling without replacement from a finite universe, J. Amer. Statist. Assoc. 47(260): 663–685.
  * Isaki and Fuller (1982) Isaki, C. T. and Fuller, W. A. (1982).  Survey design under the regression superpopulation model, J. Amer. Statist. Assoc. 77: 89–96.
  * Kang and Schafer (2007) Kang, J. D. and Schafer, J. L. (2007).  Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data, Stat. Sci. 22: 523–539.
  * Keshavan et al. (2009) Keshavan, R. H., Montanari, A. and Oh, S. (2009).  Matrix completion from noisy entries, Adv. Neural Inf. Process. Syst., Vol. 22, pp. 952–960.
  * Kim (2011) Kim, J. K. (2011).  Parametric fractional imputation for missing data analysis, Biometrika 98: 119–132.
  * Kim et al. (2006) Kim, J. K., Brick, J., Fuller, W. A. and Kalton, G. (2006).  On the bias of the multiple-imputation variance estimator in survey sampling, J. R. Stat. Soc. Ser. B. Stat. Methodol. 68: 509–521.
  * Kim and Fuller (2004) Kim, J. K. and Fuller, W. (2004).  Fractional hot deck imputation, Biometrika 91: 559–578.
  * Kim and Haziza (2014) Kim, J. K. and Haziza, D. (2014).  Doubly robust inference with missing data in survey sampling, Statist. Sinica 24(1): 375–394.
  * Kim and Yu (2011) Kim, J. K. and Yu, C. L. (2011).  A semiparametric estimation of mean functionals with nonignorable missing data, J. Amer. Statist. Assoc. 106: 157–165.
  * Koltchinskii et al. (2011) Koltchinskii, V., Lounici, K. and Tsybakov, A. B. (2011).  Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion, Ann. Statist. 39(5): 2302–2329.
  * Kott (1994) Kott, P. S. (1994).  A note on handling nonresponse in sample surveys, J. Amer. Statist. Assoc. 89(426): 693–696.
  * Kott and Chang (2010) Kott, P. S. and Chang, T. (2010).  Using calibration weighting to adjust for nonignorable unit nonresponse, J. Amer. Statist. Assoc. 105: 1265–1275.
  * Mao, Chen and Wong (2019) Mao, X., Chen, S. X. and Wong, R. K. (2019).  Matrix completion with covariate information, J. Amer. Statist. Assoc. 114(525): 198–210.
  * Mao et al. (20xx) Mao, X., Wang, Z. and Yang, S. (20xx).  Supplement to “Matrix completion for survey data prediction with multivariate missingness”.
  * Mao, Wong and Chen (2019) Mao, X., Wong, R. K. and Chen, S. X. (2019).  Matrix completion under low-rank missing mechanism, arXiv preprint arXiv:1812.07813 .
  * Mazumder et al. (2010) Mazumder, R., Hastie, T. and Tibshirani, R. (2010).  Spectral regularization algorithms for learning large incomplete matrices, J. Mach. Learn. Res. 11: 2287–2322.
  * Meng (1994) Meng, X.-L. (1994).  Multiple-imputation inferences with uncongenial sources of input, Stat. Sci. 9: 538–558.
  * Negahban and Wainwright (2012) Negahban, S. and Wainwright, M. J. (2012).  Restricted strong convexity and weighted matrix completion: optimal bounds with noise, J. Mach. Learn. Res. 13(1): 1665–1697.
  * Nielsen (2003) Nielsen, S. F. (2003).  Proper and improper multiple imputation, Int. Stat. Rev. 71: 593–607.
  * Robin et al. (2019) Robin, G., Klopp, O., Josse, J., Moulines, É. and Tibshirani, R. (2019).  Main effects and interactions in mixed and incomplete data frames, J. Amer. Statist. Assoc. (just-accepted): 1–31.
  * Robins et al. (1994) Robins, J. M., Rotnitzky, A. and Zhao, L. P. (1994).  Estimation of regression coefficients when some regressors are not always observed, J. Amer. Statist. Assoc. 89: 846–866.
  * Rubin (1976) Rubin, D. B. (1976).  Inference and missing data, Biometrika 63: 581–592.
  * Rubin (1978) Rubin, D. B. (1978).  Multiple imputations in sample surveys-a phenomenological bayesian approach to nonresponse, Proceedings of the Survey Research Methods Section of the American Statistical Association, Vol. 1, American Statistical Association, pp. 20–34.
  * Shao and Steel (1999) Shao, J. and Steel, P. (1999).  Variance estimation for survey data with composite imputation and nonnegligible sampling fractions, J. Amer. Statist. Assoc. 94: 254–265.
  * Sun and Zhang (2012) Sun, T. and Zhang, C.-H. (2012).  Calibrated elastic regularization in matrix completion, Adv. Neural Inf. Process. Syst., Vol. 25, pp. 863–871.
  * Tang et al. (2003) Tang, G., Little, R. J. and Raghunathan, T. E. (2003).  Analysis of multivariate missing data with nonignorable nonresponse, Biometrika 90: 747–764.
  * van Buuren and Groothuis-Oudshoorn (2011) van Buuren, S. and Groothuis-Oudshoorn, K. (2011).  mice: Multivariate imputation by chained equations in R, J. Stat. Softw. 45(3): 1–67.
  * van der Linden and Hambleton (2013) van der Linden, W. J. and Hambleton, R. K. (2013).  Handbook of Modern Item Response Theory, Springer, New York, NY.
  * Wang and Robins (1998) Wang, N. and Robins, J. M. (1998).  Large-sample theory for parametric multiple imputation procedures, Biometrika 85: 935–948.
  * Wang et al. (2014) Wang, S., Shao, J. and Kim, J. K. (2014).  An instrument variable approach for identification and estimation with nonignorable nonresponse, Statist. Sinica 24(3): 1097–1116.
  * Xu et al. (2013) Xu, M., Jin, R. and Zhou, Z.-H. (2013).  Speedup matrix completion with side information: application to multi-label learning, Adv. Neural Inf. Process. Syst., Vol. 26, pp. 2301–2309.
  * Yang and Kim (2016) Yang, S. and Kim, J. K. (2016).  A note on multiple imputation for method of moments estimation, Biometrika 103(1): 244–251.
  * Yang et al. (2019) Yang, S., Wang, L. and Ding, P. (2019).  Causal inference with confounders missing not at random, Biometrika, accepted .
  * Zou and Hastie (2005) Zou, H. and Hastie, T. (2005).  Regularization and variable selection via the elastic net, J. R. Stat. Soc. Ser. B. Stat. Methodol. 67(2): 301–320.