Survey Data Integration for Distribution Function Estimation

Note that if model (2.1) holds, $F_{\text{N}}(t)$ may be restated as

Now, letting

denote the finite population distribution function of $\epsilon$ at some fixed point $\frac{t-m(\boldsymbol{X}_{u};\boldsymbol{\beta})}{\nu(\boldsymbol{X}_{u})}$, we note

which further implies that

Thus, it is not egragious to replace $\mathbbm{1}\left(Y_{u}\leq t\right)$ in $F_{\text{N}}(t)$ with $G_{\text{N}}\left(\boldsymbol{X}_{u}\right)$, assuming that the working model (2.1) holds.

Given that $\boldsymbol{A}$ is subject to missingness on $Y$ only, replacing $\mathbbm{1}\left(Y_{i}\leq t\right)$ in $\widehat{F}_{\pi}(t)$ with an empirical estimate of $G_{\text{N}}(\boldsymbol{X}_{i})$ from $\boldsymbol{B}$ remains plausible if $\boldsymbol{\widehat{\beta}}\xrightarrow{P}\boldsymbol{\beta}$ (i.e., ignorability holds). To describe the procedure, let $\widehat{G}(\boldsymbol{X}_{k})$ denote the eCDF of the estimated residuals from sample $\boldsymbol{B}$, $\hat{\epsilon}_{j}=\frac{Y_{j}-m(\boldsymbol{X}_{j};\boldsymbol{\widehat{\beta}})}{\nu(\boldsymbol{X}_{j})}$, such that

Then, the residual-based estimator of $F_{\text{N}}(t)$ is defined as

Before we discuss the asymptotic properties of $\widehat{F}_{\text{R}}(t)$, we must establish a framework that allows $n_{\text{A}}$, $n_{\text{B}}$, and $N$ to be infinite. Similar to [28, 43], let $\mathscr{U}\coloneqq\set{\mathcal{U}_{\text{N}}}$ denote an increasing sequence of finite populations of size $N\to\infty$ from model (2.1); furthermore, let $\boldsymbol{A}$ and $\boldsymbol{B}$ denote samples from $\mathcal{U}_{\text{N}}$ ($N$ henceforth omitted for brevity) according to some sampling design, such that the sampling fraction $\frac{n_{\text{s}}}{N}=\frac{n_{\text{A}}+n_{\text{B}}}{N}$ converges to a limit in (0,1) as both $n_{\text{s}}$ and $N$ tend to infinity [24]. Given this asymptotic framework, let $\boldsymbol{\beta}^{*}=\text{plim}_{n_{\text{B}}\to\infty}\boldsymbol{\widehat{\beta}},$ such that, for any $\varepsilon>0$,

Letting $\xrightarrow{P}$ again denote convergence in probability, we have previously mentioned that $\boldsymbol{\widehat{\beta}}\xrightarrow{P}\boldsymbol{\beta}$ (i.e., $\boldsymbol{\beta}^{*}=\boldsymbol{\beta}$) under the ignorability condition; and since $\boldsymbol{\beta}_{\text{N}}\xrightarrow{P}\boldsymbol{\beta}$ as well, we may use Theorem 2.1.3 of [32] to show $\boldsymbol{\widehat{\beta}}\xrightarrow{P}\boldsymbol{\beta}_{\text{N}}$.

In the following, we formally define a set of conditions that govern the validity of our asymptotic results.

Assumption 1 allows $\boldsymbol{\widehat{\beta}}\xrightarrow{P}\boldsymbol{\beta}_{\text{N}}$; Assumption 2 establishes a framework to allow $n_{\text{A}},n_{\text{B}}\to\infty$ as $N\to\infty$, and limits the growth rate of both sample sizes; Assumptions 3 and 4 ensure design-based consistency under $p(\boldsymbol{A})$, the sampling design of $\boldsymbol{A}$, and allows $n_{\text{B}}$ to grow at a faster rate than $\mathbb{E}_{\mathcal{D}}\left(n_{\text{A}}\right)$; Assumption 5 ensures asymptotic normality under a general design; Assumption 6 establishes regularity conditions for $\boldsymbol{\widehat{\beta}}$; and Assumption 7 specifies a limiting smooth function for $F_{\text{N}}(t)$ [43, 24, 36]. Note that, if Assumptions 3 and 6 hold, then $\boldsymbol{\widehat{\beta}}=\boldsymbol{\beta}_{\text{N}}+o_{P}\left([\mathbb{E}_{\mathcal{D}}\left(n_{\text{A}}\right)]^{-\frac{1}{2}}\right)$.

The following theorem summarizes the asymptotic properties of the proposed residual-based CDF estimator, $\widehat{F}_{\text{R}}(t)$.

In Section 4.4, we define a replication-type variance estimator based on the bootstrap, and compare its performance with that of the variance estimator in Theorem 2.

One may consider a more direct ‘mass imputation’ approach to CDF estimation by directly substituting $Y_{i}$ for all $i\in\boldsymbol{A}$ with predicted values, $\hat{Y}_{i}\coloneqq m(\boldsymbol{X}_{i};\boldsymbol{\widehat{\beta}})$, but this approach will almost always result in higher estimation bias. Let $\widehat{F}_{\text{P}}(t)$ denote this direct plug-in CDF estimator such that

In the best-case scenario of $\boldsymbol{\widehat{\beta}}=\boldsymbol{\beta}$, we observe

and thus

which is not unbiased for $F_{\text{N}}(t)$ under model 2.1 unless $\nu(\boldsymbol{X})\epsilon_{u}$ is effectively zero for all $u\in\mathcal{U}$. From the perspective of sample $\boldsymbol{A}$, we note $\hat{Y}_{i}=Y_{i}-\varsigma_{i}$, where the unknown residual $\varsigma_{i}$ could trigger an erroneous $\mathbbm{1}(\cdot)$ output even if it is small for all $i\in\boldsymbol{A}$—that is to say, the plug-in approach is ‘all-or-nothing’, in the sense that $\mathbbm{1}(\hat{Y}_{i}\leq t)$ is either right or wrong with very little room for error. This is not the case for $\widehat{F}_{\text{R}}(t)$, though, since small $\varsigma_{i}$ implies a well-specified (and transportable) $\boldsymbol{\widehat{\beta}}$, making $\widehat{G}(\boldsymbol{X}_{i})$ a reasonable substitute for $\mathbbm{1}(Y_{i}\leq t)$ within $\widehat{F}_{\pi}(t)$.

We conclude this section with a discussion of the quantile estimation problem. From Eq. 1.2, an estimator of $t_{\text{N}}(\alpha)$ using $\widehat{F}_{\text{R}}(t)$ could be defined as

or, for $\widehat{F}_{\text{P}}(t)$, as

where $\alpha$ denotes the percentile of interest. In general, if some $\widehat{F}(t)$ is asymptotically unbiased, then its corresponding quantile, $\hat{t}(\alpha)$, will also be asymptotically unbiased [24]. And, assuming $\widehat{F}(t)$ is asymptotically normally distributed, $\gamma\%$ confidence intervals can be constructed using [44]’s formula as follows:

where $z_{\gamma/2}$ is the $(1-\gamma/2)^{th}$ percentile under the standard normal distribution and $\widehat{\text{SE}}$ is the estimated standard error of $\widehat{F}(t)$. Then, using [31, 26]’s linear approximations, we may estimate the variance of $\hat{t}(\alpha)$ using

An alternative variance estimator based on the bootstrap will be defined and empirically evaluated in Section 4.4.

One should note that quantile estimation is extremely sensitive to small sample sizes, generally, and small inclusion probabilities, specifically, as both tend to exclude many (if not most) of the quantiles in $\mathcal{U}$. In the next sections, we will evaluate the consequences of small $n_{\text{B}}$ and $n_{\text{A}}$ on the performance of both $\hat{t}_{\text{P}}(\alpha)$ and $\hat{t}_{\text{R}}(\alpha)$ using simulated data.

To thoroughly assess the performance of our proposed estimators, we considered the following four superpopulation models, $\{\xi_{i}\}^{4}_{i=1}$, which account for varying levels of complexity between $Y$ and $\boldsymbol{X}$.

• •

  Model $\xi_{1}$ [25]:

  $$Y=.3+2X_{1}+2X_{2}+\epsilon,$$

  where $X_{1},X_{2}\sim\mathrm{N}(\mu=2,\sigma=1)$ and $\epsilon\sim\mathrm{N}(\mu=0,\sigma=1)$.

• •

  Model $\xi_{2}$ [25]:

  $$Y=.3+.5X^{2}_{1}+.5X^{2}_{2}+\epsilon,$$

  where $X_{1},X_{2}\sim\mathrm{N}(\mu=2,\sigma=1)$ and $\epsilon\sim\mathrm{N}(\mu=0,\sigma=1)$.

• •

  Model $\xi_{3}$ [41, 35]:

  $$Y=-\sin{(X_{1})}+X^{2}_{2}+X_{3}-\mathrm{e}^{-X^{2}_{4}}+\epsilon,$$

  where $X_{1},\cdots,X_{4}\sim\mathrm{Uniform}\left(\min=-1,\max=1\right)$ and $\epsilon\sim\mathrm{N}\left(\mu=0,\sigma=\sqrt{.5}\right)$.

• •

  Model $\xi_{4}$ [38, 35]:

  	$\displaystyle Y$	$\displaystyle=X_{1}+.707X^{2}_{2}+2\mathbbm{1}\left(X_{3}>0\right)+.873\ln\left(|X_{1}|\right)|X_{3}|+.894X_{2}X_{4}$
  		$\displaystyle+2\mathbbm{1}\left(X_{5}>0\right)+.464\mathrm{e}^{X_{6}}+\epsilon,$

  where $X_{1},\cdots,X_{6}\sim\mathrm{N}\left(\mu=0,\sigma=1\right)$ and $\epsilon\sim\mathrm{N}\left(\mu=0,\sigma=1\right)$.

Under each of these models, we generated a finite population $\mathcal{U}$ of size $N=100,000$. After $\mathcal{U}$ was obtained, a simple random sample without replacement (SRS) of size $n_{\text{A}}=500$ was selected to serve as $\boldsymbol{A}$, the probability sample. For sample $\boldsymbol{B}$ with $\boldsymbol{n}_{\text{B}}=[500\ \ 1000\ \ 5000\ \ 10000]$, we used modifications of [25]’s stratified simple random sampling technique to generate samples under both missing at random (MAR) and missing not at random (MNAR) mechanisms. For MAR, we defined $X^{*}$ to be the covariate with the strongest correlation to $Y$, and placed all elements $u\in\mathcal{U}$ with $X_{u}^{*}\leq\mu_{X^{*}}$ in stratum I and the rest in stratum II. We then used stratified SRS to obtain $n_{\text{I}}=0.85\boldsymbol{n}_{\text{B}}$ and $n_{\text{II}}=(1-0.85)\boldsymbol{n}_{\text{B}}$ elements from each respective strata. This procedure was nearly identical to that used to generate MNAR data, with the exception of using $Y_{u}\leq\mu_{Y}$ as the strata discrimination rule. Lastly, using each of the $n_{\text{sim}}=1500$ randomly sampled $\boldsymbol{A}$/$\boldsymbol{B}$ sample pairs, we calculated $\widehat{F}_{\text{B}}(t)$ (shortened as ‘B’), $\widehat{F}_{\text{P}}(t)$ (shortened as ‘P’) and $\widehat{F}_{\text{R}}(t)$ (shortened as ‘R’), as well as their respective quantile estimators, at the 1st, 10th, 25th, 50th, 75th, 90th, and 99th percentiles. The two estimators $\widehat{F}_{\text{P}}(t)$ and $\widehat{F}_{\text{R}}(t)$ used multivariable linear regression, which was fit on sample $\boldsymbol{B}$ using R’s lm() function. Results are presented in Tables 4.1–4.1.