Matrix Completion for Survey Data Prediction with Multivariate Missingness

TL;DR

This paper introduces a novel matrix completion-based imputation method for survey data with multivariate missingness, leveraging both row and column patterns for improved accuracy, demonstrated through simulations and NHANES data application.

Contribution

It develops a new matrix completion approach that treats the data matrix holistically, incorporating covariates and a low-rank residual, with a projection strategy for parameter identification.

Findings

Doubly robust estimator achieves lower mean squared error.

Method outperforms existing imputation schemes in simulations.

Applied successfully to NHANES data.

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.

Tables

Table 1. 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 𝑛 500 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

Table 2. 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 𝐿 500 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

Table 3. 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

Equations

$${\mathbf{Y}}_N={\mathbf{A}}_N+{\mathbf{ϵ}}_N,$$

$${\mathbf{A}}_N={\mathbf{X}}_N{\mathbf{\beta }}^{\ast }+{\mathbf{B}}_N^{\ast },$$

$$p_{ij}=p_{ij}\left({\mathbf{x}}_i\right)=\frac{\exp \left\{\left(1,{\mathbf{x}}_i^{\mathrm{T}}\right){\mathbf{\gamma }}_{.j}\right\}}{1+\exp \left\{\left(1,{\mathbf{x}}_i^{\mathrm{T}}\right){\mathbf{\gamma }}_{.j}\right\}},$$

$$R(\mathbf{\beta },\mathbf{B})=\frac{1}{NL}E{\Vert {\mathbf{R}}_N\circ {\mathbf{P}}_N^{†}\circ {\mathbf{Y}}_N-{\mathbf{X}}_N\mathbf{\beta }-\mathbf{B}\Vert }_F^2,$$

$${\hat{\theta }}_j=\frac{1}{N}\sum _{i=1}^N\frac{I_i}{{\pi }_i}{y}_{ij}.$$

$$R^{\ast }(\mathbf{\beta },{\mathbf{B}}_n)$$

$$=$$

$${\mathbf{D}}_n^{-1/2}\left({\mathbf{R}}_n\circ {\mathbf{P}}_n^{†}\circ {\mathbf{Y}}_n\right)-{\mathbf{D}}_n^{-1/2}{\mathbf{X}}_n\mathbf{\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).$$

$${\mathbf{\beta }}^{\ast \prime }={\mathbf{\beta }}^{\ast }+({\mathbf{X}}_n^{\mathrm{T}}{\mathbf{D}}_n^{-1}{\mathbf{X}}_n{)}^{-1}{\mathbf{X}}_n^{\mathrm{T}}{\mathbf{D}}_n^{-1}{\mathbf{B}}_n^{\ast }\text{and}{\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 }),$$

$$R^{\ast }(\mathbf{\beta },{\mathbf{B}}_n)=R^{\ast }({\mathbf{\beta }}^{\prime },{\mathbf{B}}_n^{\prime })=$$

$${\Vert {\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^{†}\circ {\mathbf{Y}}_n\right)\right\}-{\mathbf{B}}_n^{\prime }\Vert }_F^2].$$

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

$${\hat{R}}^{\ast }({\mathbf{\beta }}^{\prime },{\mathbf{B}}_n^{\prime })=$$

$${\Vert {\mathcal{P}}_{{\mathbf{D}}_n^{-1/2}{\mathbf{X}}_n}^{\perp }\left\{{\mathbf{D}}_n^{-1/2}\left({\mathbf{R}}_n\circ {\hat{\mathbf{P}}}_n^{†}\circ {\mathbf{Y}}_n\right)\right\}-{\mathbf{B}}_n^{\prime }\Vert }_F^2],$$

$$({\hat{\mathbf{\beta }}}^{\prime },{\hat{\mathbf{B}}}_n^{\prime })=\underset{{\mathbf{\beta }}^{\prime },{\mathbf{B}}_n^{\prime }\in \mathcal{N}({\mathbf{D}}_n^{-1/2}{\mathbf{X}}_n)}{\mathrm{argmin}}{\hat{R}}^{\ast }({\mathbf{\beta }}^{\prime },{\mathbf{B}}_n^{\prime })+{\tau }_1{\Vert {\mathbf{\beta }}^{\prime }\Vert }_F^2+{\tau }_2\left\{\alpha {\Vert {\mathbf{B}}_n^{\prime }\Vert }_{\ast }+(1-\alpha ){\Vert {\mathbf{B}}_n^{\prime }\Vert }_F^2\right\},$$

$${\hat{\mathbf{\beta }}}^{\prime }={\left({\mathbf{X}}_n^{\mathrm{T}}{\mathbf{D}}_n^{-1}{\mathbf{X}}_n+NL{\tau }_1\mathbf{I}\right)}^{-1}{\mathbf{X}}_n^{\mathrm{T}}{\mathbf{D}}_n^{-1}\left({\mathbf{R}}_n\circ {\hat{\mathbf{P}}}_n^{†}\circ {\mathbf{Y}}_n\right).$$

$${\mathcal{T}}_c\left(\mathbf{M}\right)=\mathbf{U}\mathrm{diag}(\{{\left({\sigma }_i-c\right)}_{+}\}){\mathbf{V}}^{⊺}$$

$${\hat{\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 {\hat{\mathbf{P}}}_n^{†}\circ {\mathbf{Y}}_n\right)\right\}\right].$$

$${\hat{\mathbf{A}}}_n={\mathbf{X}}_n{\hat{\mathbf{\beta }}}^{\prime }+{\mathbf{D}}_n^{1/2}{\hat{\mathbf{B}}}_n^{\prime }.$$

$${\hat{\theta }}_{j,DR}=\frac{1}{N}\sum _{i=1}^N\frac{I_i}{{\pi }_i}\left\{\frac{r_{ij}(y_{ij}-{\hat{a}}_{ij})}{{\hat{p}}_{ij}}+{\hat{a}}_{ij}\right\},$$

$${\hat{\theta }}_{j,DR}=\frac{1}{N}\sum _{i=1}^N\frac{I_i}{{\pi }_i}\left\{\frac{r_{ij}(y_{ij}-{\hat{a}}_{ij})}{p_{ij}}+{\hat{a}}_{ij}\right\}+o_P(1)$$

$$\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\},$$

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

$$\frac{1}{mL}{\Vert {\hat{\mathbf{\beta }}}^{\prime }-{\mathbf{\beta }}^{\ast \prime }\Vert }_F^2\le C_1{r}_{{\mathbf{B}}_N}{L}^{-1}\log \left(n\right)p_{\min }^{-1}\text{and}$$

$$\frac{1}{nL}{\Vert {\hat{\mathbf{B}}}_n^{\prime }-{\mathbf{B}}_n^{\ast \prime }\Vert }_F^2\le C_2{r}_{{\mathbf{B}}_N}N{n}^{-1}{L}^{-1}\log \left(n\right)p_{\min }^{-1}.$$

$${\hat{\mathbf{A}}}_n={\mathbf{X}}_n{\hat{\mathbf{\beta }}}^{\prime }+{\mathbf{D}}_n^{1/2}{\hat{\mathbf{B}}}_n^{\prime },$$

$$\frac{1}{nL}{\Vert {\hat{\mathbf{A}}}_n-{\mathbf{A}}_n\Vert }_F^2\le C_3{r}_{{\mathbf{B}}_N}{L}^{-1}\log \left(n\right)p_{\min }^{-1}.$$

$$\frac{1}{L}\sum _{j=1}^L({\hat{\theta }}_{j,DR}-{\theta }_j{)}^2=O_p\{r_{{\mathbf{B}}_N}{L}^{-1}\log \left(n\right)\}.$$

$${\hat{\theta }}_{j,IPM}=N^{-1}\sum _{i=1}^n{r}_{ij}{\hat{p}}_{ij}^{-1}{y}_{ij}.$$

$$y_{ij}={\varphi }_{0j}+{\mathbf{x}}_i^{\mathrm{T}}{\mathbf{\varphi }}_{1j},$$

$$({\hat{\varphi }}_{0j},{\hat{\mathbf{\varphi }}}_{1j})=\underset{({\varphi }_{0j},{\mathbf{\varphi }}_{1j})}{\mathrm{argmin}}\sum _{i=1}^n\frac{r_{ij}}{{\pi }_i{\hat{p}}_{ij}}(y_{ij}-{\varphi }_{0j}-{\mathbf{x}}_i^{\mathrm{T}}{\mathbf{\varphi }}_{1j}{)}^2.$$

$$\text{mbox}{Bias}_j={\hat{\theta }}_j^{(m)}-{\theta }_j\text{and}\text{mbox}{SE}_j=\frac{1}{1000}\sum _{m=1}^{1000}({\hat{\theta }}_j^{(m)}-{\hat{\theta }}_j{)}^2,$$

Full text

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$,

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

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}$,

[TABLE]

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):

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math] 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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE] 4. IV

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

[TABLE]

and the parameters in (13) are estimated by

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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.

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.

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).