Matrix completion under complex survey sampling

Abstract

Multivariate nonresponse is often encountered in complex survey sampling, and simply ignoring it leads to erroneous inference. In this paper, we propose a new matrix completion method for complex survey sampling. Different from existing works either conducting row-wise or column-wise imputation, the data matrix is treated as a whole which allows for exploiting both row and column patterns simultaneously. A column-space-decomposition model is adopted incorporating a low-rank structured matrix for the finite population with easy-to-obtain demographic information as covariates. Besides, we propose a computationally efficient projection strategy to identify the model parameters under complex survey sampling. Then, an augmented inverse probability weighting estimator is used to estimate the parameter of interest, and the corresponding asymptotic upper bound of the estimation error is derived. Simulation studies show that the proposed estimator has a smaller mean squared error than other competitors, and the corresponding variance estimator performs well. The proposed method is applied to assess the health status of the U.S. population.

Appendices

Appendix

A Technical conditions

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

• C1 (a) The random errors $\{\epsilon _{ij}:i=1,\ldots ,N;j=1,\ldots ,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$, ${\max }_{i,j}E|\epsilon _{ij}|^{l}\le \frac{1}{2}l!c_{\sigma }^{2}\eta ^{l-2}$ for any positive integer $l\ge 2$.

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

• C3 The population design matrix ${\varvec{X}}_N$ is of size $N\times d$ such that $N>d$. Moreover, there exists a positive constant $a_{x}$ such that $\Vert {\varvec{X}}_N\Vert _{\infty }\le a_{x}$ and ${\varvec{X}}_N^{\mathrm{T}}{\varvec{D}}_N{\varvec{X}}_N$ is invertible, where ${\varvec{D}}_N$ is a diagonal matrix with $\pi _i$ as its (i, i)th entry. Furthermore, there exists a symmetric matrix ${\varvec{S}}_{{\varvec{X}}}$ with $\sigma _{\min }({\varvec{S}}_{{\varvec{X}}})\asymp 1\asymp \Vert {\varvec{S}}_{{\varvec{X}}}\Vert$ such that $n_0^{-1}{\varvec{X}}_N^{\mathrm{T}}{\varvec{D}}_N{\varvec{X}}_N\rightarrow {\varvec{S}}_{{\varvec{X}}}$ as $N\rightarrow \infty$, where $n_0=\sum _{i=1}^N\pi _i$ is the expected sample.

• C4 There exists a positive constant a such that $\max \{\Vert {\varvec{X}}_N\varvec{\beta }^{*}\Vert _{\infty },\Vert {\varvec{A}}_N\Vert _{\infty }\}\le a$.

• C5 The indicators of observed entries $\{r_{ij}:i=1,\ldots ,N;j=1,\ldots ,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 ${\varvec{X}}_N$. Furthermore, for $i=1,\dots ,N$ and $j=1,\dots ,L$, $\Pr (r_{ij}=1 | {\varvec{x}}_{i}, y_{ij}) = \Pr (r_{ij}=1 | {\varvec{x}}_{i})$ follows the logistic regression model (7).

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

• C7 The sampling design satisfies that $N^{-1}\sum _{i=1}^Ny_i\pi _i^{-1} = O_p(n^{-1/2})$ if $N^{-1}\sum _{i=1}^Ny_i^{2}$ is asymptotically bounded.

Condition C1(a) is a common regularity condition for the measurement errors in $\varvec{\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). 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}{\varvec{X}}_N^{\mathrm{T}}{\varvec{D}}_N{\varvec{X}}_N$ is fixed at $d\times d$, it is mild to assume ${\varvec{X}}_N^{\mathrm{T}}{\varvec{D}}_N{\varvec{X}}_N$ to be invertible, and there exists a symmetric matrix ${\varvec{S}}_{{\varvec{X}}}$ as the limit of $n_0^{-1}{\varvec{X}}_N^{\mathrm{T}}{\varvec{D}}_N{\varvec{X}}_N$. Please notice that we do not assume randomness for generating ${\varvec{X}}_N$, and it is a common assumption for design-based framework. 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}{\varvec{X}}_n^{\mathrm{T}}{\varvec{X}}_n$ is also ${\varvec{S}}_{{\varvec{X}}}$ under regularity conditions. The order of $\sigma _{\min }({\varvec{S}}_{{\varvec{X}}})$ and $\Vert {\varvec{S}}_{{\varvec{X}}}\Vert$ equals to 1 is due to $\Vert {\varvec{X}}_N\Vert _{\infty }<\infty$. Condition C4 is also standard in the matrix completion literature (Cai and Zhou, 2016; Koltchinskii et al., 2011; Negahban and Wainwright, 2012). 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 (7). 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\le n$ is quite mild. Condition C7 is a mild restriction on the estimator for the population mean, and it can be satisfied under general sampling designs. To get general results, we do not make any assumptions for the asymptotic relationship between the population size N and the sample size n; see Theorem 1 for details. We can make further assumptions for the sample sizes to guarantee certain convergence properties; see the discussion of Theorem 3.

B Lemmas

Under the logistic model (7), together with the results in Mao et al. (2019) and Sweeting (1980), it can been shown that for all $t>d+3$, there exist some positive constants $C_g$, C and $C_{d}$ such that $\Pr \{\sum _{ij}(1/\widehat{p}_{ij}-1/p_{ij})^{2}\ge C_g p_{\min }^{-3}t\}\le C_{d}t\exp \{-t/2\}+L\max _{j}\sup _{t}|\Pr \{\sum _{i}(1/\widehat{p}_{ij}-1/p_{ij})^{2}\ge t\}-\Pr (\chi ^{2}_{d+1}\ge C_g p_{\min }^{-3}t)|$. Then, $\max _{j}\sup _{t}|\Pr \{\sum _{i}(1/\widehat{p}_{ij}-1/p_{ij})^{2}\ge t\}-\Pr (\chi ^{2}_{d+1}\ge C_g p_{\min }^{-3}t)|\le L^{-2}$ and $C_{d}t\exp \{-t/2\}$ is a function independent of n and L such that $\lim _{t\rightarrow \infty }t\exp \{-t/2\}= 0$.

Write ${\varvec{J}}_{ij}={\varvec{e}}_{i}(n_1){\varvec{e}}^\intercal _{j}(n_2)$, where ${\varvec{e}}_{i}(n)\in {\mathbb {R}}^n$ is the standard basis vector with the i-th element being 1 and the rest being 0. Now we present several lemmas.

Lemma 1

Under Conditions C2 and C3 and Poisson sampling, we have

$$\begin{aligned} n^{-1}{\varvec{X}}_n^{\mathrm{T}}{\varvec{X}}_n={\varvec{S}}_{{\varvec{X}}}+ o_p(1). \end{aligned}$$

(15)

Proof of Lemma 1

Denote ${\varvec{e}}_i$ to be a column vector of length d with jth element being 1 and others being 0. Recall that $n_0 = \sum _{k=1}^N\pi _k$ is the expected sample size. For $i=1,\ldots ,d$ and $j=1,\ldots ,d$, consider

$$\begin{aligned} E( n_0^{-1}{\varvec{e}}_i^{\mathrm{T}}{\varvec{X}}_n^{\mathrm{T}}{\varvec{X}}_n{\varvec{e}}_j)&= n_0^{-1}E\left( \sum _{k=1}^NI_kx_{ki}x_{kj}\right) = n_0^{-1}\sum _{i=1}^Nx_{ki}\pi _kx_{kj} \nonumber \\&= n_0^{-1}{\varvec{e}}_i^{\mathrm{T}}{\varvec{X}}_N^{\mathrm{T}}{\varvec{D}}_N{\varvec{X}}_N{\varvec{e}}_j, \end{aligned}$$

(16)

where the expectation is take with respect to the sampling design, and ${\varvec{x}}_i$ is the ith row of ${\varvec{X}}_N$.

Under Poisson sampling, we have

$$\begin{aligned} \mathrm {var}( n_0^{-1}{\varvec{e}}_i^{\mathrm{T}}{\varvec{X}}_n^{\mathrm{T}}{\varvec{X}}_n{\varvec{e}}_j)&= n_0^{-2}\sum _{k=1}^N\pi _k(1-\pi _k)x_{k,i}x_{k,j}\nonumber \\ &< n_0^{-2}\sum _{k=1}^N\pi _kx_{k,i}x_{k,j}\nonumber \\ &= n_0^{-1}\left( n_0^{-1}{\varvec{e}}_i^{\mathrm{T}}{\varvec{X}}_N^{\mathrm{T}}{\varvec{D}}_N{\varvec{X}}_N{\varvec{e}}_j\right) . \end{aligned}$$

(17)

By a similar argument in Wang et al. (2019), we can show that $n/n_0\rightarrow 1$ in probability as $N\rightarrow \infty$ under Condition C2. By Condition C3, (16) and (17), we have proved Lemma 1. $\square$

Lemma 2

Let $\Psi ^{(1)}=\sum _{ij}r_{ij}\epsilon _{ij}{\varvec{J}}_{ij}/(nL\widehat{p}_{ij}\pi _{ij}^{1/2})$. Under Conditions C1, C5 and C6 and Poisson sampling, for some positive constants $C_1$, $c_{\sigma }$, $\eta$, $\delta _{\sigma }$ and all $t>d+3$, we have

$$\begin{aligned}&\left\Vert \Psi ^{(1)}\right\Vert \\&\quad \le C_1\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\} \end{aligned}$$

holds with probability at least $1-1/(n+L)-C_{d}t\exp \{-t/2\}-1/L-12c_{\sigma }^{2}\eta ^{2}\log ^{-\delta _{\sigma }}(n)$.

Lemma 3

Let $\Psi ^{(2)}=\sum _{ij}a_{ij}(r_{ij}/p_{ij}-1){\varvec{J}}_{ij}/(nL\pi _{ij}^{1/2})$. Under Conditions C4–C6 and Poisson sampling, for some positive constants $C_2$, we have

$$\begin{aligned} \left\Vert \Psi ^{(2)}\right\Vert \le C_2N^{1/2}n^{-1}L^{-1}\log ^{1/2}\left( n\right) p_{\min }^{-1/2} \end{aligned}$$

holds with probability at least $1-1/(n+L)$.

Lemma 4

Let $\Psi ^{(3)}=\sum _{ij}a_{ij}(r_{ij}/{\widehat{p}}_{ij}-r_{ij}/p_{ij}){\varvec{J}}_{ij}/(nL\pi _{ij}^{1/2})$. Under Conditions C4 and C6 and Poisson sampling, for some positive constants $C_3$, $\delta _{\sigma }$ and all $t>d+3$, we have

$$\begin{aligned} \left\Vert \Psi ^{(3)}\right\Vert \le C_3N^{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} \end{aligned}$$

holds with probability at least $1-C_{d}t\exp \{-t/2\}-1/L$.

It is easy to show Lemma 2–4 by the proof of Lemma S4.1–S4.3 in the supplementary material of Mao et al. (2019).

C Proofs

1.1 C.1 Proof of Theorem 1

With the definition of $\Delta (\delta _{\sigma },t)$ in (12), under Conditions C1–C6 and Poisson sampling, together with Lemmas 2–4, We have for a positive constant $C_0$,

$$\begin{aligned} \left\Vert \Psi ^{(1)}\right\Vert +\Vert \Psi ^{(2)}\Vert +\left\Vert \Psi ^{(3)}\right\Vert \le C_{0} \Delta (\delta _{\sigma },t), \end{aligned}$$

with probability at least $1-2/n-2C_{d}t\exp \{-t/2\}-2/L-12c_{\sigma }^{2}\eta ^{2}\log ^{-\delta _{\sigma }}(n)$.

Let ${\varvec{X}}_n^{\prime }={\varvec{D}}_n^{-1/2}{\varvec{X}}_n$ and $\eta _{n,L}(\delta _{\sigma },t) = 4/(n+L)+ 4C_{d}t\exp \{-t/2\}+4/L+C\log ^{-\delta _{\sigma }}(n)$ for a positive constant C. By choosing t as (13), $\tau _1\asymp N^{-1}nL^{-1}\log ^{-1/2}(n)\Delta (\delta _{\sigma })$ and $\tau _2\asymp \eta _{g}^{-1/2}N^{-1}n^{1/4}L^{-1/4}\log ^{1/2}(L)\log ^{\delta _{\sigma }/3}(n)$, where $\Delta (\delta _{\sigma })=N^{1/2}n^{-1}L^{-1}\log ^{1/2}(n)p_{\min }^{-1/2}$ and $1-\alpha \asymp (nL)^{-1}$ in (8) for any $\delta _{\sigma }>0$, together with Condition C2 and Poisson sampling, it follows the same proof with the proof of Corollary 1 in Mao et al. (2019) that, for some constants $C_1$ and $C_2$, with probability at least $1-\eta _{n,L}(\delta _{\sigma },t)$,

$$\begin{aligned} \text {both}\quad \frac{1}{nL}\left\Vert {\varvec{X}}_n{\widehat{\varvec{\beta }}}^{\prime }-{\varvec{X}}_n\varvec{\beta }^{*\prime }\right\Vert _F^{2} \quad \text {and}\quad \frac{1}{nL}\left\Vert {\widehat{{\varvec{B}}}}_n^{\prime }-{\varvec{B}}_n^{*\prime }\right\Vert _F^{2}\le C_1r_{{\varvec{B}}_N} Nn^{-1}L^{-1}\log \left( n\right) p_{\min }^{-1}. \end{aligned}$$

Thus it is easy to obtain that

$$\begin{aligned} \frac{1}{mL}\left\Vert {\widehat{\varvec{\beta }}}^{\prime }-\varvec{\beta }^{*\prime }\right\Vert _F^{2}\le C_2r_{{\varvec{B}}_N} L^{-1}\log \left( n\right) p_{\min }^{-1}, \end{aligned}$$

under Condition C3. $\square$

1.2 C.2 Proof of Theorem 2

Due to the observations that

$$\begin{aligned} \left\Vert {\widehat{{\varvec{A}}}}_n-{\varvec{A}}_n\right\Vert _{F}^2\le \left\Vert {\varvec{X}}_n{\widehat{\varvec{\beta }}}^{\prime }-{\varvec{X}}_n\varvec{\beta }^{*\prime }\right\Vert _F^{2}+\left\Vert {\varvec{D}}_n^{-1/2}{\widehat{{\varvec{B}}}}_n^{\prime }-{\varvec{D}}_n^{-1/2}{\varvec{B}}_n^{*\prime }\right\Vert _F^{2}, \end{aligned}$$

together with Theorem 1, it is easy to obtain the result under Condition C2 and Poisson sampling.$\square$

1.3 C.3 Proof of Theorem 3

Denote

$$\begin{aligned}&{\tilde{\theta }}_{j,{\textit{AIPW}}} = N^{-1}\sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij} - {a}_{ij})}{{p}_{ij}} + {a}_{ij}\right\} , \end{aligned}$$

(18)

$$\begin{aligned}&{\theta }^\dag _{j,{\textit{AIPW}}} = N^{-1}\sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij} - \widehat{a}_{ij})}{{p}_{ij}} + \widehat{a}_{ij}\right\} \end{aligned}$$

(19)

for $j=1,\ldots ,L$. The difference between $\widehat{\theta }_{j,{\textit{AIPW}}}$ in (10) and ${\tilde{\theta }}_{j,{\textit{AIPW}}}$ in (18) is that we use estimators $\widehat{N}$, $\widehat{p}_{ij}$ and $\widehat{a}_{ij}$ for $\widehat{\theta }_{j,{\textit{AIPW}}}$ but use true values N, $p_{ij}$ and $a_{ij}$ for ${\tilde{\theta }}_{j,{\textit{AIPW}}}$. The difference between $\widehat{\theta }_{j,{\textit{AIPW}}}$ and ${\theta }^\dag _{j,{\textit{AIPW}}}$ in (19) is that we use $\widehat{N}$ for $\widehat{\theta }_{j,{\textit{AIPW}}}$, but use N for ${\theta }^\dag _{j,{\textit{AIPW}}}$.

First, we prove

$$\begin{aligned} {\tilde{\theta }}_{j,{\textit{AIPW}}} - \theta _j = O_p(n^{-1/2}). \end{aligned}$$

(20)

Consider

$$\begin{aligned} E( {\tilde{\theta }}_{j,{\textit{AIPW}}}) = E\{E({\tilde{\theta }}_{j,{\textit{AIPW}}}\mid \{I_i\})\} = N^{-1}\sum _{i=1}^N\frac{E(I_i)}{\pi _i}(y_{ij}) = \theta _j, \end{aligned}$$

(21)

where the first equality holds due to $E(r_{ij}) =p_{ij}$. Next, we derive the variance of ${\tilde{\theta }}_{j,{\textit{AIPW}}}$. Specifically, we have

$$\begin{aligned} \mathrm {var}({\tilde{\theta }}_{j,{\textit{AIPW}}})&= \frac{1}{N^2}E\left( \mathrm {var}\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}\right\} \mid S\right] \right) \nonumber \\&\quad + \frac{1}{N^2}\mathrm {var}\left( E\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}\right\} \mid S\right] \right) = V_{1,j} + V_{2,j}, \end{aligned}$$

(22)

where $S=\{I_i:i=1,\ldots ,N\}$.

Because $E(r_{ij}) = p_{ij}$, we have

$$\begin{aligned} \mathrm {var}\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}\right\} \mid S\right] = \sum _{i=1}^N\frac{I_i(1-p_{ij})}{\pi _i^2p_{ij}}(y_{ij}-a_{ij})^2. \end{aligned}$$

Thus, we have

$$\begin{aligned} V_{1,j} = \frac{1}{N^2}\sum _{i=1}^N\frac{1-p_{ij}}{\pi _{i}}(y_{ij}-a_{ij})^2=O_p(n^{-1}), \end{aligned}$$

(23)

where the last equality holds by Conditions C1, C2 and the strong law of large numbers (Athreya and Lahiri, 2006). Notice that

$$\begin{aligned} E\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}\right\} \mid S\right] =\sum _{i=1}^N\frac{I_iy_{ij}}{\pi _i}. \end{aligned}$$

By the models (1)–(2) and Condition C4, we can show that $N^{-1}\sum _{i=1}^Ny_i^2$ is asymptotically bounded. Thus, by Condition C7, we have

$$\begin{aligned} V_{2,j}=O_p(n^{-1}) \end{aligned}$$

(24)

By (21)–(24), we have shown (20).

Next, we show that

$$\begin{aligned} \frac{1}{L}\sum _{j=1}^L({\theta }^\dag _{j,{\textit{AIPW}}} - {\tilde{\theta }}_{ j,{\textit{AIPW}}})^2=O_p\{r_{{\varvec{B}}_N}L^{-1}\log \left( n\right) \}. \end{aligned}$$

(25)

Consider

$$\begin{aligned} {\theta }^\dag _{j,{\textit{AIPW}}} - {\tilde{\theta }}_{ j,{\textit{AIPW}}} = \frac{1}{N}\sum _{i=1}^N\frac{I_i}{\pi _{i}}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})(p_{ij}-\widehat{p}_{ij})}{p_{ij}\widehat{p}_{ij}} + \frac{(r_{ij}-\widehat{p}_{ij})(a_{ij}-\widehat{a}_{ij})}{\widehat{p}_{ij}}\right\} . \end{aligned}$$

(26)

Consider

$$\begin{aligned} E\left\{ \frac{1}{N}\sum _{i\in S}\pi _{i}^{-1}\frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} \right\}&= \frac{1}{N}\sum _{i=1}^N(y_{ij}-a_{ij})=O_p(N^{-1/2}), \end{aligned}$$

(27)

$$\begin{aligned} \mathrm {var}\left\{ \frac{1}{N}\sum _{i\in S}\pi _{i}^{-1}\frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} \right\}&= O_p(n^{-1}), \end{aligned}$$

(28)

where the asymptotic order in (28) holds due to Condition C7 and $N^{-1}\sum _{i=1}^N(y_{ij}-a_{ij})^2$ is asymptotically bounded in probability since $\{\epsilon _{ij}:j=1,\ldots ,N\}$ are independent and their variances are uniformly bounded. By (27) and (28), we have

$$\begin{aligned} \frac{1}{N}\sum _{i=1}^N\frac{I_i}{\pi _{i}}\frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} = O_p(n^{-1/2} ). \end{aligned}$$

(29)

Because the response model (7) for $p_{ij}$ is assumed to be correctly specified, and $p_{ij}$ is bounded away from 0 by Condition C6, we have $\widehat{p}_{ij}^{-1}(p_{ij}-\widehat{p}_{ij}) = O_p(1)$ uniformly for $i=1,\ldots ,N$. Thus, by (29), we have

$$\begin{aligned} N^{-1}\sum _{i=1}^N\frac{I_i}{\pi _{i}}\frac{r_{ij}(y_{ij}-a_{ij})(p_{ij}-\widehat{p}_{ij})}{p_{ij}\widehat{p}_{ij}} =O_p(n^{-1/2}). \end{aligned}$$

(30)

By (26) and (30), we have

$$\begin{aligned} {\theta }^\dag _{j,{\textit{AIPW}}} - {\tilde{\theta }}_{ j,{\textit{AIPW}}} = O_p(n^{-1/2}) + \frac{1}{N}\sum _{i=1}^N\frac{I_i}{\pi _i}\frac{(r_{ij}-\widehat{p}_{ij})(a_{ij}-\widehat{a}_{ij})}{\widehat{p}_{ij}}. \end{aligned}$$

(31)

Since $\widehat{p}_{ij}=p_{ij}+o_p(1)$, we have

$$\begin{aligned} \frac{r_{ij}-\widehat{p}_{ij}}{\widehat{p}_{ij}}=\{1+o_p(1)\}\frac{r_{ij}-p_{ij}}{p_{ij}}. \end{aligned}$$

(32)

Since $p_{ij}\ge p_{\mathrm {min}}>0$ by Condition C6, we have

$$\begin{aligned} \frac{r_{ij}-\widehat{p}_{ij}}{\widehat{p}_{ij}}= O_p(1) \end{aligned}$$

(33)

uniformly for $i=1,\ldots ,N$.

Thus, by Condition C2, (31) and (33), we have

$$\begin{aligned} \frac{1}{L}\sum _{j=1}^L({\theta }^\dag _{j,{\textit{AIPW}}} - {\tilde{\theta }}_{ j,{\textit{AIPW}}})^2&\le O_p(n^{-1})+ \frac{C_6 O_p(1)}{Ln^2}\sum _{j=1}^L\left\{ \sum _{i=1}^n(a_{ij} - \widehat{a}_{ij}) \right\} ^2,\nonumber \\ &\le O_p(n^{-1})+ \frac{ O_p(1)}{Ln}\sum _{j=1}^L\sum _{i=1}^n(a_{ij} - \widehat{a}_{ij})^2\nonumber \\ &=O_p(n^{-1})+ \frac{O_p(1)}{nL}\left\Vert \widehat{A}_n-A_n\right\Vert _F^2, \end{aligned}$$

(34)

where $\pi _i\ge C_6^{-1}nN^{-1}$ for $C_6>0$ by Condition C2, and we have assumed that the first n subjects are sampled. By (20), (34) and Theorem 2 and the fact that $L\le n$, we have proved (25).

By Condition C4 and the fact that $E(\epsilon ^2_{ij})<\sigma _0^2$ uniformly, $\theta _j$ is uniformly bounded for $j=1,\ldots ,L$ in probability. Thus, by (20) and (25), we conclude that

$$\begin{aligned} \frac{1}{L}\sum _{j=1}^L({\theta }^\dag _{j,{\textit{AIPW}}})^2 = O_p\{r_{{\varvec{B}}_N}L^{-1}\log \left( n\right) \} \end{aligned}$$

(35)

By Condition C7, we conclude that $\widehat{N}N^{-1}=1+O_p(n^{-1/2})$. Consider

$$\begin{aligned} \frac{1}{L}\sum _{j=1}^L(\widehat{\theta }_{j,{\textit{AIPW}}} - {\theta }^\dag _{ j,{\textit{AIPW}}})^2 = \frac{O_p(n^{-1})}{L}\sum _{j=1}^L({\theta }^\dag _{ j,{\textit{AIPW}}})^2 = o_p\{r_{{\varvec{B}}_N}L^{-1}\log \left( n\right) \}, \end{aligned}$$

(36)

where the first equality holds since $\widehat{N}N^{-1}=1+O_p(n^{-1/2})$ uniformly for $j=1,\ldots ,L$, and the second equality holds by (35). Thus, by (34) and (36), we have proved Theorem 3. $\square$

D Plug-in variance estimators

When deriving the plug-in variance estimator, we ignore the variability for estimating $\widehat{a}_{ij}$. First, we consider the plug-in variance estimator under Poisson sampling. For $j=1,\ldots ,L$, let

$$\begin{aligned} g_j(\theta ) = \frac{1}{N}\sum _{i=1}^{N}\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-{a}_{ij})}{{p}_{ij}}+{a}_{ij}-\theta \right\} \end{aligned}$$

be the estimating function for the AIPW estimator with $\widehat{a}_{ij}$ replaced by ${a}_{ij}$. Let ${\tilde{\theta }}_{j,{\textit{AIPW}}}$ solves $g_j(\theta )=0$, and we use a variance estimator of ${\tilde{\theta }}_{j,{\textit{AIPW}}}$ to approximate that of $\widehat{\theta }_{j,{\textit{AIPW}}}$.

It can be shown that ${\tilde{\theta }}_{j,{\textit{AIPW}}}-\theta _j=O_p(n^{-1/2})$, so we have

$$\begin{aligned} 0 = g_j({\tilde{\theta }}_{j,{\textit{AIPW}}}) = g_j(\theta _j) + g_j'(\theta _j) ({\tilde{\theta }}_{j,{\textit{AIPW}}}-\theta _j) + o_p(n^{-1/2}), \end{aligned}$$

(37)

where $g_j'(\theta ) = -N^{-1}\sum _{i=1}^NI_i\pi _i^{-1}$ is the derivative of $g_j(\theta )$. By a similar argument for (27)–(28), we can show that $g_j'(\theta _j)\rightarrow -1$ in probability. Besides, by (37), we have

$$\begin{aligned} ({\tilde{\theta }}_{j,{\textit{AIPW}}}-\theta _j) = -\{g_j'(\theta _j)\}^{-1}g_j(\theta _j) + o_p(n^{-1/2}) = g_j(\theta _j)+ o_p(n^{-1/2}). \end{aligned}$$

Thus, the variance of ${\tilde{\theta }}_{j,{\textit{AIPW}}}$ can be estimated by the one of $g_j(\theta _j)$.

Consider

$$\begin{aligned} \mathrm {var}\{g_j(\theta _j)\}&= N^{-2}E\left( \mathrm {var}\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}-\theta _j\right\} \mid S\right] \right) \nonumber \\&\quad + N^{-2}\mathrm {var}\left( E\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}-\theta _j\right\} \mid S\right] \right) \nonumber \\&= V_{1,j} + V_{2,j}. \end{aligned}$$

(38)

Since $E(r_{ij}) = p_{ij}$, we have

$$\begin{aligned} \mathrm {var}\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{m_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}-\theta _j\right\} \mid S\right] = \sum _{i=1}^N\frac{I_i(1-p_{ij})}{\pi _i^2p_{ij}}(y_{ij}-a_{ij})^2. \end{aligned}$$

Thus, we have

$$\begin{aligned} V_{1,j} = N^{-2}\sum _{i=1}^N\frac{1-p_{ij}}{\pi _ip_{ij}}(y_{ij}-a_{ij})^2, \end{aligned}$$

(39)

and it can be estimated by

$$\begin{aligned} \widehat{V}_{1,j} = N^{-2}\sum _{i=1}^N\frac{r_{ij}(1-p_{ij})}{\pi _i^2p_{ij}^2}(y_{ij}-a_{ij})^2. \end{aligned}$$

(40)

Notice that

$$\begin{aligned} E\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}-\theta _j\right\} \mid S\right] =\sum _{i=1}^N\frac{I_i}{\pi _i}(y_{ij}-\theta _j). \end{aligned}$$

Under Poisson sampling,

$$\begin{aligned} \mathrm {var}\left( E\left[ \sum _{i=1}^N\frac{I_i}{\pi _i}\left\{ \frac{r_{ij}(y_{ij}-a_{ij})}{p_{ij}} + a_{ij}-\theta _j\right\} \mid S\right] \right) = \sum _{i=1}^N\frac{1-\pi _i}{\pi _i}(y_{ij}-\theta _j)^2. \end{aligned}$$

Thus,

$$\begin{aligned} V_{2,j} = N^{-2} \sum _{i=1}^N\frac{1-\pi _i}{\pi _i}(y_{ij}-\theta _j)^2, \end{aligned}$$

(41)

and it can be estimated by

$$\begin{aligned} \widehat{V}_{2,j} = N^{-2}\sum _{i=1}^N\frac{I_ir_{ij}(1-\pi _i)}{p_{ij}\pi _i^2}y_{ij}^2. \end{aligned}$$

(42)

By (38)–(42) and plugging in $\widehat{a}_{ij}$ and $\widehat{\theta }_{j,{\textit{AIPW}}}$ for $a_{ij}$ and $\theta _j$, respectively, the plug-in variance estimator for ${\widehat{\theta }}_{j,{\textit{AIPW}}}$ is

$$\begin{aligned} \widehat{V}_{j,poi}&= N^{-2}\sum _{i=1}^N\frac{I_ir_{ij}(1-p_{ij})}{\pi _i^2p_{ij}^2}(y_{ij}-\widehat{a}_{ij})^2\\&\quad + N^{-2}\sum _{i=1}^N\frac{I_ir_{ij}(1-\pi _i)}{p_{ij}\pi _i^2}(y_{ij}-\widehat{\theta }_{j,{\textit{AIPW}}})^2 \end{aligned}$$

under Poisson sampling.

Use a similar argument, we can show that the plug-in variance estimator is

$$\begin{aligned} {\hat{V}}_{ j,srs}&= n^{-2}\sum _{i=1}^n\frac{I_ir_{ij}(1-{\hat{p}}_{ij})}{{\hat{p}}_{ij}^2}(y_{ij}-{\widehat{a}}_{ij})^2+ n^{-1}(1-nN^{-1})\\&\quad \left[ n^{-1}\sum _{i=1}^N\frac{I_ir_{ij}(y_{ij}-\widehat{\theta }_{j,{\textit{AIPW}}})^2}{{\hat{p}}_{ij}} - \left( n^{-1}\sum _{i=1}^N\frac{I_ir_{ij}(y_{ij}-\widehat{\theta }_{j,{\textit{AIPW}}})}{{\hat{p}}_{ij}}\right) ^2 \right] \end{aligned}$$

under simple random sampling, and the one is

$$\begin{aligned} {\hat{V}}_{ j,pps}&= N^{-2}\sum _{i\in S}\frac{r_{ij}(1-{\hat{p}}_{ij})}{(nq_i)^2{\hat{p}}_{ij}^2}(y_{ij}-{\widehat{a}}_{ij})^2 + \{n(n-1)\}^{-1}\\&\quad \left\{ \sum _{i\in S}\frac{r_{ij}(y_{ij}-\widehat{\theta }_{j,{\textit{AIPW}}})^2}{{\hat{p}}_{ij}q_i^2} - n^{-1}\left( \sum _{i\in S}\frac{r_{ij}(y_{ij}-\widehat{\theta }_{j,{\textit{AIPW}}})}{{\hat{p}}_{ij}q_i}\right) ^2\right\} \end{aligned}$$

under probability-proportional-to-size sampling, where S is the index set of the sample, and $q_i$ is the selection probability of the ith element.

E Balanced repeated replication method

Consider a stratified multi-stage sampling design with two clusters selected per stratum for the first stage. Denote $w_{hik}$ to be the survey weight associated with $y_{hik}$, the kth sample element in the ith cluster of the hth stratum. The basic idea of the modified balanced repeated replication is to use the same estimation method based on the “reconstruct” the survey weight, that is,

$$\begin{aligned} w_{hik}^{(r)} = w_{hik} (1+\epsilon \delta _{rh}), \end{aligned}$$

where $\epsilon \in (0,1)$ is a predefined constant, and $\delta _{rh}=1$ or $\delta _{rh} = -1$ for the rth repetition. A set of R repetitions is said to be balanced if $\sum _{r=1}^R\delta _{rh}\delta _{rh'}=0$ for $h\ne h'$. The $R\times H$ matrix $(\delta _{rh})_{R\times H}$ can be obtained from a Hadamard matrix. Please check Rao and Shao (1999) for details about the modified balanced repeated replication method. In the simulation study and real data application, we choose $\epsilon = 1/2$ and $R =H$.