Doubly Robust Estimation of the Finite Population Distribution Function Using Nonprobability Samples

1. Introduction

Probability sampling, first introduced in Neyman’s seminal work [1], has long been regarded as the gold standard for finite population inference in survey statistics [2]. In contrast, nonprobability samples were historically employed mainly in observational studies or as supplementary sources for probability sampling. In recent years, however, probability sampling has encountered growing challenges, including rising survey costs, incomplete sampling frames, declining response rates, and increasingly strict data privacy regulations [3,4]. These difficulties have renewed interest in nonprobability samples as a practical alternative for population inference [5,6].

A defining feature of nonprobability samples is their unknown selection mechanism, which, if ignored, can lead to substantial selection bias [7,8,9]. Moreover, the theoretical foundations for inference with nonprobability samples remain underdeveloped, and standardized criteria for evaluating the reliability of resulting estimates are largely lacking. A prominent strategy to address these limitations is data integration, which leverages information from high-quality probability samples to adjust for the selection bias inherent in nonprobability samples. While early research on data integration focused primarily on combining two probability samples [10,11], more recent work has expanded the scope to include methods that incorporate diverse nonprobability sources, such as online panels and opt-in datasets [12,13].

Within the data integration framework, several approaches have been developed to mitigate selection bias. For example, Elliott and Valliant [13] proposed an approach that models the selection mechanism of the nonprobability sample and applies the inverse of estimated propensity scores as weights. Alternatively, Kim et al. [14] introduced a mass imputation method, which fits an outcome regression model using study variables observed in the nonprobability sample and predicts corresponding values for units in the probability sample. A central limitation of both approaches, however, is that their validity depends entirely on correct model specification.

To overcome these limitations, doubly robust estimation methods have attracted considerable attention, as they yield asymptotically unbiased estimates if at least one of the two underlying models is correctly specified [6,15,16,17]. Existing applications of doubly robust estimation, however, have largely concentrated on measures of central tendency, such as the population mean, with relatively little attention to distributional estimands—such as finite distribution function and quantiles—that are essential for characterizing the overall shape of a population distribution. While the mean is intuitive and easily interpretable, it is highly sensitive to extreme values and may be inadequate when information about specific regions of the distribution is crucial, for instance in studies of income inequality or health outcomes. By contrast, finite distribution function and quantiles provide information across the entire distribution, enabling more comprehensive analyses [18,19].

In this study, we extend the doubly robust estimator for the population mean proposed by Chen et al. [16] to the estimation of the finite distribution function within the data integration framework. Specifically, we develop three estimators based on auxiliary variables observed in both probability and nonprobability samples: an inverse probability weighted estimator, a regression estimator, and a doubly robust estimator. The proposed doubly robust estimator attains asymptotic unbiasedness for the finite distribution function, provided that either the propensity score model or the outcome regression model is correctly specified.

The remainder of this paper is organized as follows. Section 2 presents the basic setup and assumptions. Section 3 defines the three estimators of the finite distribution function and analyzes their theoretical properties. Section 4 applies these estimators to the estimation of population quantiles and describes the construction of Woodruff confidence intervals. Section 5 evaluates the performance of the proposed methods through simulation studies. Finally, Section 6 concludes with the implications of the study and directions for future research.

2. Basic Setup

Consider a finite population $U=\{1,2,\dots ,N\}$ of size N. Each unit $i\in U$ is associated with auxiliary variables ${\mathbf{x}}_i$ and a study variable $y_i$. The parameter of interest in this study is the finite population distribution function, defined as

$$F_y\left(t\right)={\displaystyle \frac{1}{N}}\sum _{i\in U}I(y_i\le t),-\infty <t<\infty ,$$

(1)

where $I\left(A\right)$ is the indicator function that equals 1 if the event A is true and 0 otherwise.

Suppose two samples are drawn from the population. The first is a nonprobability sample $S_A$ of size $n_A$, in which both the auxiliary variables ${\mathbf{x}}_i$ and the study variable $y_i$ are observed for each unit in $i\in S_A$. The second is a probability sample $S_B$ of size $n_B$, in which the auxiliary variables ${\mathbf{x}}_i$ are observed together with their inclusion probabilities ${\pi }_i^B=P(i\in S_B)$ under a given sampling design p. In this setup, the nonprobability sample is not representative of the population, whereas the probability sample does not contain observations on the study variable. An integrated approach is therefore required for population inference. This approach relies on common auxiliary variables observed in both samples, which serve as a bridge linking the study variable with the design information. Such a data integration framework is well-established and has been widely applied in prior studies [6,12,13,16,17,20,21]. Following Chen et al. [16], we adopt this framework and assume that the two samples are drawn independently, which simplifies the analysis and enhances theoretical validity.

To utilize nonprobability samples for population inference, we assume that they are drawn under an implicit probabilistic selection mechanism, even though their inclusion probabilities are unknown. This corresponds to the concept of an unknown probability sample discussed by Wu [22]. Under this assumption, the issue of undercoverage is excluded from the scope of this study.

To formalize the selection mechanism, define the indicator variable $R_i$ as

$$R_i=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}i\in S_A,\hfill \\ 0,\hfill & \mathrm{if}i\notin S_A.\hfill \end{array}\right.$$

The conditional expectation of $R_i$ is

$${\pi }_i^A={\mathrm{E}}_{\delta }[R_i\mid {\mathbf{x}}_i,y_i]=P(R_i=1\mid {\mathbf{x}}_i,y_i),$$

(2)

where the subscript $\delta$ denotes the model for the selection mechanism. This probability ${\pi }_i^A$ is commonly referred to as the propensity score in observational studies [23] and as the participation probability in survey sampling [6,24].

We adopt the concept of an unknown probability sample and impose the following assumptions on the propensity score, as introduced in Chen et al. [16].

A1  

Given the auxiliary variables ${\mathbf{x}}_i$, the study variable $y_i$ and the selection indicator $R_i$ are conditionally independent.

A2  

For all units $i\in U$, ${\pi }_i^A>0$.

A3  

Given ${\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N$, the selection indicators $R_1,R_2,\dots ,R_N$ are independent.

Assumption A1 and A2 together constitute the strong ignorability condition [23]. Under this condition, ${\pi }_i^A$ in Equation (2) depends only on the auxiliary variables ${\mathbf{x}}_i$.

3. Estimators of the Finite Distribution Function

In this section, we introduce three estimators for the finite distribution function in Equation (1) under the data integration framework: the inverse probability weighted (IPW) estimator, the regression (REG) estimator, and the doubly robust (DR) estimator. Their statistical properties are examined within a joint randomization framework consisting of the following three components [17]:

(i)

$\delta$: selection mechanism for the nonprobability sample

(ii)

p: the probability sampling design

(iii)

M: the regression model for the study variable y

Specifically, the IPW estimator is analyzed under the $\delta p$-framework, the REG estimator under the $Mp$-framework, and the DR estimator under either the $\delta p$- or the $Mp$-framework, without specification of which one. All these frameworks incorporate the probability sampling design p.

Following Chen et al. [16], we adopt conditions C1, C4, C5, and C6, which are redefined as Regularity Conditions B1–B4 in this paper, modified by substituting $I(y_i\le t)$ for $y_i$ and $G(t-{\mathbf{x}}_i^{\top }\beta )$ for $m({\mathbf{x}}_i,\beta )$. To ensure the validity of the Taylor expansion, we also impose the following additional condition

B5  

The distribution function $G\left(t\right)$ of the error term is twice differentiable for all t.

For asymptotic analysis, we consider a sequence of finite populations indexed by $\nu$, denoted as $U_{\nu }:\nu =1,2,\dots$ with corresponding samples $S_{A,\nu }$ and $S_{B,\nu }$. As $\nu \to \infty$, the population size $N_{\nu }$ and the sample sizes $n_{A,\nu },n_{B,\nu }$ all diverge to infinity. For simplicity, the subscript $\nu$ is omitted hereafter, and asymptotics are expressed in terms of $N\to \infty$.

3.1. Inverse Probability Weighted Estimator

The IPW method is widely used to correct for selection bias in nonprobability samples, where the inclusion probability ${\pi }_i^A$ is assumed to be a function of auxiliary variables ${\mathbf{x}}_i$ and unknown parameters $\theta$ of a participation model

$${\pi }_i^A=\pi ({\mathbf{x}}_i,\theta )$$

Because auxiliary information for the entire population is often unavailable, Chen et al. [16] proposed a pseudo-likelihood approach that incorporates the design weights from a probability sample. The participation model parameters are first estimated by $\widehat{\theta }$, which are then used to compute estimated inclusion probabilities ${\widehat{\pi }}_i^A=\pi ({\mathbf{x}}_i,\widehat{\theta })$. The corresponding pseudo-weight is ${\widehat{d}}_i^A=1/{\widehat{\pi }}_i^A$, which is used to construct a Hájek-type estimator of the finite distribution function. This quasi-randomization approach enables valid inference from nonprobability samples [13,25].

For a fixed point t, the IPW estimator of the finite distribution function is defined as

$${\widehat{F}}_{\mathrm{IPW}}\left(t\right)={\displaystyle \frac{1}{{\widehat{N}}_A}}\sum _{i\in S_A}{\widehat{d}}_i^{A}I(y_i\le t),$$

(3)

where ${\widehat{N}}_A={\sum }_{i\in S_A}{\widehat{d}}_i^A$.

Result 1.  

UnderA1–A3andB1–B5, and if the propensity score model is correctly specified as a logistic regression model, ${\widehat{F}}_{\mathrm{IPW}}\left(t\right)$ is an asymptotically $\delta p$-unbiased estimator of $F_y\left(t\right)$ for fixed point t.

Proof. 

The asymptotic property of the IPW estimator for the population mean ${\mu }_y$, denoted as ${\widehat{\mu }}_y={\widehat{N}}_A^{-1}{\sum }_{i\in S_A}{y}_i/{\widehat{\pi }}_i^A$, is given by Chen et al. [16] as ${\widehat{\mu }}_y-{\mu }_y=O_p\left(n_A^{-1/2}\right)$. The finite distribution function estimator ${\widehat{F}}_{\mathrm{IPW}}\left(t\right)$ can be viewed as a plug-in estimator with $y_i$ replaced by $I(y_i\le t)$. Therefore,

$${\widehat{F}}_{\mathrm{IPW}}\left(t\right)-F_y\left(t\right)=O_p\left(n_A^{-1/2}\right),-\infty <t<\infty$$

Moreover, since both ${\widehat{F}}_{\mathrm{IPW}}\left(t\right)$ and $F_y\left(t\right)$ lie in $[0,1]$, we have $|{\widehat{F}}_{\mathrm{IPW}}\left(t\right)-F_y\left(t\right)|\le 1$. Hence,

$${\mathrm{E}}_{\delta p}\left[{\widehat{F}}_{\mathrm{IPW}}\left(t\right)-F_y\left(t\right)\right]=O\left(n_A^{-1/2}\right)$$

Thus, under the $\delta p$-randomization framework, ${\widehat{F}}_{\mathrm{IPW}}\left(t\right)$ is an asymptotically $\delta p$-unbiased estimator of the finite distribution function $F_y\left(t\right)$.    □

The IPW estimator is asymptotically unbiased under the assumption of strong ignorability, that is, when the selection mechanism is fully explained by the auxiliary variables, a condition also referred to as missing at random (MAR). However, if the propensity score model is misspecified, asymptotic unbiasedness may fail. Moreover, even under correct specification, extreme propensity score values (close to 0 or 1) may yield highly unstable inverse probability weights, inflating the variance of the estimator [26].

3.2. Regression Estimator

The REG estimator fits an outcome regression of $y_i$ on ${\mathbf{x}}_i$ using the nonprobability sample $S_A$, and then applies the fitted model to units in the probability sample $S_B$ to estimate the finite distribution function. Because it imputes the study variable for all units in $S_B$, this procedure is commonly referred to as mass imputation. This approach is based on a superpopulation model, in which the finite population $U=\{1,2,\dots ,N\}$ is assumed to be a random sample generated from the following outcome regression model

$$y_i={\mathbf{x}}_i^{\top }\beta +{ϵ}_i,i=1,2,\dots ,N,$$

(4)

where the error terms ${ϵ}_i$ follow a normal distribution with ${\mathrm{E}}_M\left[{ϵ}_i\right]=0$ and ${\mathrm{Var}}_M\left[{ϵ}_i\right]={\sigma }^2$, and are assumed to be independent. The subscript M indicates that the corresponding expectation or variance is taken under the model (4). By Assumption A1, we have ${\mathrm{E}}_M[y_i\mid {\mathbf{x}}_i,R_i]={\mathrm{E}}_M[y_i\mid {\mathbf{x}}_i]$ and ${\mathrm{Var}}_M[y_i\mid {\mathbf{x}}_i,R_i]={\mathrm{Var}}_M[y_i\mid {\mathbf{x}}_i]$, which ensures that the model is valid for the nonprobability sample as well.

A simple model-based approach to estimating the finite distribution function replaces $I(y_i\le t)$ with the plug-in indicator $I({\mathbf{x}}_i^{\top }\widehat{\beta }\le t)$. However, since ${\mathrm{E}}_M\left[I({\mathbf{x}}_i^{\top }\widehat{\beta }\le t)\right]\ne P(y_i<t)$ in general, this naive substitution can lead to bias [18]. To address this, Chambers and Dunstan [27] proposed a residuals-based estimator. Let $G(·)$ denote the distribution function of ${ϵ}_i$. The residual-based estimator is then defined as

$${\widehat{G}}_i{\left(t\right)}_{\mathrm{REG}}={\displaystyle \frac{1}{n_A}}\sum _{j\in S_A}I(y_j-{\mathbf{x}}_j^{\top }\widehat{\beta }\le t-{\mathbf{x}}_i^{\top }\widehat{\beta }),$$

(5)

where $\widehat{\beta }$ is the ordinary least squares estimator of $\beta$.

Based on Equation (5), the regression estimator of $F_y\left(t\right)$ at a fixed point t is given by

$${\widehat{F}}_{\mathrm{REG}}\left(t\right)={\displaystyle \frac{1}{{\widehat{N}}_B}}\sum _{i\in S_B}{d}_i^B{\widehat{G}}_i{\left(t\right)}_{\mathrm{REG}},$$

(6)

where ${\widehat{N}}_B={\sum }_{i\in S_B}{d}_i^B$.

Result 2. 

UnderA1–A3andB1–B5, ${\widehat{F}}_{\mathrm{REG}}\left(t\right)$ is an asymptotically $Mp$-unbiased estimator of $F_y\left(t\right)$ for fixed point t.

Proof. 

Following the lines of Chambers et al. [28], it can be shown that the result holds under conditions B1–B5

$$E_M\left[{\widehat{G}}_i{\left(t\right)}_{\mathrm{REG}}\right]=G(t-{\mathbf{x}}_i^{\top }\beta )+O\left(n_A^{-1}\right).$$

Using this result, the expectation of the bias can be evaluated as

$$\begin{array}{cc}\hfill \mathrm{E}\left[{\widehat{F}}_{\mathrm{REG}}\left(t\right)-F_y\left(t\right)\right]& ={\mathrm{E}}_p\left[{\mathrm{E}}_M\left[{\widehat{F}}_{\mathrm{REG}}\left(t\right)-F_y\left(t\right)\right]\right]\hfill \\ \hfill & ={\mathrm{E}}_p\left[{\displaystyle \frac{1}{{\widehat{N}}_B}}\sum _{i\in S_B}{d}_i^B\left\{G(t-{\mathbf{x}}_i^{\top }\beta )+O\left(n_A^{-1}\right)\right\}-{\displaystyle \frac{1}{N}}\sum _{i=1}^{N}G(t-{\mathbf{x}}_i^{\top }\beta )\right]\hfill \\ \hfill & \approx {\displaystyle \frac{1}{N}}\sum _{i=1}^{N}G(t-{\mathbf{x}}_i^{\top }\beta )-{\displaystyle \frac{1}{N}}\sum _{i=1}^{N}G(t-{\mathbf{x}}_i^{\top }\beta )+O\left(n_A^{-1}\right)\hfill \\ \hfill & =O\left(n_A^{-1}\right)\hfill \end{array}$$

Therefore, under the $Mp$-randomization framework, ${\widehat{F}}_{\mathrm{REG}}\left(t\right)$ is an asymptotically $Mp$-unbiased estimator for the finite distribution function $F_y\left(t\right)$.    □

When the outcome model is correctly specified, the REG estimator is highly efficient and supports broader use of nonprobability samples. However, if the regression model fails to capture the true distribution, bias may arise and the method becomes sensitive to misspecification. To mitigate this limitation, the next section introduces the DR estimator, which remains asymptotically unbiased provided that either the propensity score model or the outcome regression model is correctly specified.

3.3. Doubly Robust Estimator

The asymptotic unbiasedness of the IPW estimator in Equation (3) and REG estimator in Equation (6) hinges on correct specification of their respective working models. In practice, however, such correctness is difficult to guarantee, motivating procedures that are robust to model misspecification. The DR estimator was introduced to address this issue and has been regarded as a successful approach since Robins et al. [29].

To construct a DR estimator for the finite distribution function, we require an analogue of ${\widehat{G}}_i{\left(t\right)}_{\mathrm{REG}}$ in Equation (5) that estimates the error distribution $G(·)$ and remains valid under the joint randomization. Because ${\widehat{G}}_i{\left(t\right)}_{\mathrm{REG}}$ is derived under the $Mp$-framework, it cannot be directly applied when the selection mechanism $\delta$ is operative (i.e., under the $\delta p$-framework). Accordingly, we extend the method of Rao et al. [30] and propose a new estimator of the error distribution that is valid under such joint randomization.

$${\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}={\displaystyle \frac{1}{{\widehat{N}}_A}}\sum _{j\in S_A}{\widehat{d}}_j^{A}I(y_j-{\mathbf{x}}_j^{\top }\widehat{\beta }\le t-{\mathbf{x}}_i^{\top }\widehat{\beta })$$

(7)

Based on ${\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}$ defined in Equation (7), the DR estimator of the finite distribution function $F_y\left(t\right)$ at a fixed point t is then given by

$${\widehat{F}}_{\mathrm{DR}}\left(t\right)={\displaystyle \frac{1}{{\widehat{N}}_A}}\sum _{i\in S_A}{\widehat{d}}_i^A\left\{I(y_i\le t)-{\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}\right\}+{\displaystyle \frac{1}{{\widehat{N}}_B}}\sum _{i\in S_B}{d}_i^B{\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}$$

(8)

Result 3. 

Under regularity conditionsA1–A3andB1–B5, and if at least one of the propensity score model or the outcome regression model is correctly specified, $F_{\mathrm{DR}}\left(t\right)$ is an asymptotically unbiased estimator of $F_y\left(t\right)$ at a fixed point t under the $\delta p$- or $Mp$-framework.

Proof. 

(i)

When the propensity score model is correctly specified

The doubly robust estimator can be rewritten as

$${\widehat{F}}_{\mathrm{DR}}\left(t\right)={\widehat{F}}_{\mathrm{IPW}}\left(t\right)-{\displaystyle \frac{1}{{\widehat{N}}_A}}\sum _{i\in S_A}{\widehat{d}}_i^A{\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}+{\displaystyle \frac{1}{{\widehat{N}}_B}}\sum _{i\in S_B}{d}_i^B{\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}$$

The second and third terms on the right-hand side are Hájek estimators of $N^{-1}{\sum }_{i=1}^N{\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}$ based on the nonprobability sample $S_A$ and the probability sample $S_B$, respectively, and hence cancel out asymptotically. Given the asymptotic $\delta p$-unbiasedness of ${\widehat{F}}_{\mathrm{IPW}}\left(t\right)$, ${\widehat{F}}_{\mathrm{DR}}\left(t\right)$ is also asymptotically $\delta p$-unbiased.

(ii)

When the outcome regression model is correctly specified

Similarly to the proof for the REG estimator, we have

$${\mathrm{E}}_M\left[{\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}\right]=G(t-{\mathbf{x}}_i^{\top }\beta )+O\left(n_A^{-1}\right)$$

Using this, the expected bias is

$$\begin{array}{cc}\hfill & \mathrm{E}\left[{\widehat{F}}_{\mathrm{DR}}\left(t\right)-F_y\left(t\right)\right]\hfill \\ & ={\mathrm{E}}_p\left[{\displaystyle \frac{1}{{\widehat{N}}_A}}\sum _{i\in S_A}{\widehat{d}}_i^A{\mathrm{E}}_M\left[I(y_i\le t)-{\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}\right]+{\displaystyle \frac{1}{{\widehat{N}}_B}}\sum _{i\in S_B}{d}_i^B{\mathrm{E}}_M\left[{\widehat{G}}_i{\left(t\right)}_{\mathrm{DR}}\right]-{\displaystyle \frac{1}{N}}\sum _{i\in U}{\mathrm{E}}_M\left[I(y_i\le t)\right]\right]\hfill \\ \hfill & ={\mathrm{E}}_p\left[{\displaystyle \frac{1}{{\widehat{N}}_B}}\sum _{i\in S_B}{d}_i^{B}G(t-{\mathbf{x}}_i^{\top }\beta )-{\displaystyle \frac{1}{N}}\sum _{i=1}^{N}G(t-{\mathbf{x}}_i^{\top }\beta )+O\left(n_A^{-1}\right)\right]\hfill \\ \hfill & \approx {\displaystyle \frac{1}{N}}\sum _{i=1}^{N}G(t-{\mathbf{x}}_i^{\top }\beta )-{\displaystyle \frac{1}{N}}\sum _{i=1}^{N}G(t-{\mathbf{x}}_i^{\top }\beta )+O\left(n_A^{-1}\right)\hfill \\ \hfill & =O\left(n_A^{-1}\right)\hfill \end{array}$$

Therefore, under the $Mp$-framework, ${\widehat{F}}_{\mathrm{DR}}\left(t\right)$ is an asymptotically unbiased estimator of the finite distribution function $F_y\left(t\right)$.    □

The asymptotic unbiasedness of the DR estimator requires the estimated coefficients to satisfy probability-limit conditions; specifically, for the propensity score parameters $\widehat{\theta }$ and the outcome regression parameters $\widehat{\beta }$, there exist fixed vectors ${\theta }^{*}$ and ${\beta }^{*}$ such that $plim\widehat{\theta }={\theta }^{*}$ and $plim\widehat{\beta }={\beta }^{*}$. If the propensity score model is correctly specified, then ${\theta }^{*}=\theta$, and if the outcome regression model is correctly specified, then ${\beta }^{*}=\beta$. Under misspecification, by contrast, these probability limits need not coincide with the true parameters, and the limiting value itself does not have a meaningful interpretation.

4. Quantile Estimation

An important application of the finite distribution function estimators is the estimation of population quantiles, defined as

$${\xi }_q=inf\{t;F_y\left(t\right)\ge q\}$$

Quantiles provide informative summaries of distributional features such as central tendency, spread, and asymmetry, and they are useful for assessing the presence of outliers. Because estimators of the finite distribution function are typically step functions, linear interpolation is employed to obtain a unique estimate of the qth quantile  [30,31,32]. The quantile estimator ${\widehat{\xi }}_q$ is expressed as

$${\widehat{\xi }}_q=a+{\displaystyle \frac{q-\widehat{F}\left(a\right)}{\widehat{F}\left(b\right)-\widehat{F}\left(a\right)}}(b-a),$$

where $a=max\{t;\widehat{F}\left(t\right)\le q\}$ and $b=min\{t;\widehat{F}\left(t\right)\ge q\}$.

A widely used method for constructing a confidence interval (CI) for a quantile estimator was proposed by Woodruff [31]. The key idea is to first obtain a CI for the estimated finite distribution function and then invert this interval to derive a CI for the quantile. The resulting $100(1-\alpha )\%$ CI is given by

$$\begin{array}{c}\hfill {\widehat{\xi }}_q^L=inf\left\{t;\widehat{F}\left(t\right)\ge q-z_{1-\alpha /2}\sqrt{\widehat{\mathrm{V}}\left[\widehat{F}\left({\widehat{\xi }}_q\right)\right]}\right\},\\ \hfill {\widehat{\xi }}_q^U=inf\left\{t;\widehat{F}\left(t\right)\ge q+z_{1-\alpha /2}\sqrt{\widehat{\mathrm{V}}\left[\widehat{F}\left({\widehat{\xi }}_q\right)\right]}\right\},\end{array}$$

where $z_{1-\alpha /2}$ is the $(1-\alpha /2)$ quantile of the standard normal distribution and $\widehat{\mathrm{V}}\left[\widehat{F}\left({\widehat{\xi }}_q\right)\right]$ denotes the estimated variance of $\widehat{F}\left(t\right)$ evaluated at ${\widehat{\xi }}_q$. Sitter and Wu [33] provided empirical evidence that the Woodruff method attains approximately correct coverage even for extreme quantiles (large or small q).

5. Simulation Studies

To evaluate the performance of the proposed finite distribution function estimators, ${\widehat{F}}_{\mathrm{IPW}}\left(t\right)$, ${\widehat{F}}_{\mathrm{REG}}\left(t\right)$, and ${\widehat{F}}_{\mathrm{DR}}\left(t\right)$, we conducted simulation studies based on two populations: A synthetic finite population from Chen et al. [16], and the 2023 Korean Survey of Household Finances and Living Conditions.

The variances of the finite distribution function estimators were obtained via a bootstrap procedure following Chen et al. [16]:

1.

From the nonprobability sample $S_A$ and the probability sample $S_B$, draw bootstrap samples $S_A^{\left(j\right)}$ and $S_B^{\left(j\right)}$ of sizes $n_A$ and $n_B$, respectively, by simple random sampling with replacement, for $J=1,000$ replicates.

2.

For each bootstrap replicate, compute ${\widehat{F}}^{\left(j\right)}\left(t\right)$.

3.

Using ${\left\{F^{\left(j\right)}\left(t\right)\right\}}_{j=1}^J$ calculate the bootstrap variance estimator $v_{\mathrm{BT}}$.

Performance was then assessed over $R=3,000$ simulation replications using percentage relative bias (%RB) and relative root mean squared error (RRMSE), where

$$\%\mathrm{RB}={\displaystyle \frac{1}{R}}\sum _{r=1}^R{\displaystyle \frac{{\widehat{\theta }}^{\left(r\right)}-\theta }{\theta }}\times 100,\mathrm{RRMSE}={\displaystyle \frac{1}{\theta }}\sqrt{{\displaystyle \frac{1}{R}}{\sum }_{r=1}^R{\left({\widehat{\theta }}^{\left(r\right)}-\theta \right)}^2},$$

with ${\widehat{\theta }}^{\left(r\right)}$ denoting the estimate from replication r and $\theta$ the target parameter. For the finite distribution function, the bootstrap variance, and quantiles, the corresponding choices were

• The finite distribution function: ${\widehat{\theta }}^{\left(r\right)}={\widehat{F}}^{\left(r\right)}\left(t\right),\theta =F_y\left(t\right)$
• Bootstrap variance: ${\widehat{\theta }}^{\left(r\right)}=v_{\mathrm{BT}}^{\left(r\right)},\theta =V$
• Quantile: ${\widehat{\theta }}^{\left(r\right)}={\widehat{\xi }}_q^{\left(r\right)},\theta ={\xi }_q$

where V denotes the simulation-based variance of $\widehat{F}\left(t\right)$ computed from 10,000 replications.

The coverage probability of the confidence interval based on the bootstrap variance ($\%{\mathrm{CP}}_v$) was evaluated as

$$\%{\mathrm{CP}}_v={\displaystyle \frac{1}{R}}\sum _{r=1}^{R}I\left(|{\widehat{F}}^{\left(r\right)}\left(t\right)-F_y\left(t\right)|\le z_{1-\alpha /2}\sqrt{v_{\mathrm{BT}}^{\left(r\right)}}\right).$$

The performance of the Woodruff confidence interval was assessed by its coverage probability ($\%\mathrm{CP}\xi$), lower error rate (%L), and upper error rate (%U):

$$\begin{array}{cc}\hfill \%{\mathrm{CP}}_{\xi }& ={\displaystyle \frac{1}{R}}\sum _{r=1}^{R}I\left({\xi }_q^{L\left(r\right)}<{\xi }_q<{\xi }_q^{U\left(r\right)}\right),\hfill \\ \hfill \%\mathrm{L}& ={\displaystyle \frac{1}{R}}\sum _{r=1}^{R}I\left({\xi }_q<{\xi }_q^{L\left(r\right)}\right),\hfill \\ \hfill \%\mathrm{U}& ={\displaystyle \frac{1}{R}}\sum _{r=1}^{R}I\left({\xi }_q^{U\left(r\right)}<{\xi }_q\right),\hfill \end{array}$$

where ${\xi }_q^{L\left(r\right)}$ and ${\xi }_q^{U\left(r\right)}$ denote, respectively, the lower and upper Woodruff CI bounds for the qth quantile in replication r.

5.1. Study 1

Following the simulation design of Chen et al. [16], we generated a finite population of size $N=20,000$. The study variable y and auxiliary variables $\mathbf{x}$ are generated from

$$y_i=2+x_{1i}+x_{2i}+x_{3i}+x_{4i}+\sigma {ϵ}_i,i=1,2,\dots ,N,$$

where $(x_{1i},x_{2i},x_{3i},x_{4i})$ follow the design in Chen et al. [16], and the error terms ${ϵ}_i\sim N(0,1)$. The parameter $\sigma$ is chosen such that the correlation coefficient $\rho$ between y and the linear predictor ${\mathbf{x}}^{\top }\beta$ equals 0.5.

We consider four model specification scenarios

• TT: Both $\delta$ and M are correctly specified.
• TF: $\delta$ is correctly specified, but M is misspecified, with $x_{4i}$ omitted from the model.
• FT: M is correctly specified, but $\delta$ is misspecified, with $x_{4i}$ omitted from the model.
• FF: Both models are misspecified, with $x_{4i}$ omitted in each model.

The analysis uses a nonprobability sample $S_A$ of size $n_A=500$ and a probability sample $S_B$ of size $n_B=1,000$. Table 1 reports %RB and RRMSE for the proposed finite distribution function estimators. Under TT, all estimators exhibit low bias and error, indicating stable performance. Under TF and FT, the DR estimator attains lower bias and error than the alternatives, highlighting the advantages of the doubly robust property. By contrast, under FF, performance deteriorates substantially for all estimators.

Table 2 compares the bootstrap variance estimators in terms of %RB and $\%{\mathrm{CP}}_v$. Under TT, all variance estimators perform satisfactorily. Under TF and FT, despite model misspecification, the variance estimator associated with the DR method retains low bias and an $\%{\mathrm{CP}}_v$ close to 95%, indicating stable reliability and accuracy. Conversely, under FF, coverage performance deteriorates markedly across all methods.

Table 3 summarizes the results for the quantile estimators. Mirroring the findings for the finite distribution function estimators, all methods perform well under the TT scenario. Under TF and FT, the DR-based quantiles remain stable, confirming the robustness of the doubly robust approach. By contrast, under FF, overall estimation accuracy deteriorates.

Table 4 reports the Woodruff CI results for the quantile estimators, including $\%\mathrm{CP}\xi$, %L, and %U.Consistent with previous findings, all methods perform well under the TT scenario. Under TF and FT, the DR-based intervals maintain $\%\mathrm{CP}\xi$ close to the nominal 95% with balanced tail errors, indicating high reliability. By contrast, under FF, coverage deteriorates substantially across methods. $\%\mathrm{CP}\xi$ falls below the nominal level and both tail error rates increase, signaling degraded interval performance.

5.2. Study 2

In the second simulation study, we treat the 2023 Korean Survey of Household Finances and Living Conditions (SHFLC; $(N=16,730)$) as the finite population and repeatedly draw subsamples from it. Table 5 summarizes the key variables used in the experiment and their definitions.

The nonprobability sample $S_A$ was generated to mimic structures commonly observed in practice. The propensity score model was specified as a logistic regression,

$$log\left\{{\displaystyle \frac{{\pi }_i^A}{1-{\pi }_i^A}}\right\}={\zeta }_0+{\zeta }_1\mathtt{EDU}+{\zeta }_2\mathtt{SNG}+{\zeta }_3\mathtt{APT}+{\zeta }_4\mathtt{DEBT},i=1,\dots ,N,$$

where ${\zeta }_0$ was chosen so that ${\sum }_{i=1}^N{\pi }_i^A=n_A$. Under this design, households with higher educational attainment of the household head, non-single households, apartment residents, and households without debt were more likely to be included in $S_A$. The nonprobability sample $S_A$ was then selected by Poisson sampling with inclusion probabilities ${\pi }_i^A$.

The probability sample $S_B$ was stratified into nine strata defined by GEO, HOME, and SIZE. A mixed allocation scheme—combining Neyman and proportional allocation—was used to determine stratum specific sample sizes, followed by simple random sampling without replacement within each stratum. The sample sizes were set to $n_A=500$ and $n_B=1000$.

The study variable of interest was current income (INCOME). Because the true outcome model was unknown, we included EXP1 and EXP2- the covariates with comparatively strong explanatory power- as regressors in the working model. This setup allows us to assess the impact of model misspecification on estimation performance and to isolate efficiency gains attributable to the DR estimator. We consider two scenarios regarding the propensity score model:

• A: correctly specified propensity score model.
• B: misspecified propensity score model (excluding SNG and DEBT).

Table 6 reports the results for the distribution–function estimators. Overall, the REG estimator performs reasonably well, although its bias and error are somewhat larger at lower quantiles than at middle and upper quantiles, likely reflecting the limited explanatory power of the auxiliary variables and the possible over-representation of high-income households. Under Scenario A, the IPW estimator and the DR estimator both exhibit low bias and error, confirming the effectiveness of propensity score adjustment. Under Scenario B, REG estimator is the most stable, while the DR estimator inherits some bias from the misspecified IPW component and thus loses efficiency. In summary, when the propensity score model is correctly specified, the IPW estimator, the REG estimator, and the DR estimator all yield stable results. However, when the propensity-score model is misspecified, only the REG estimator and the DR estimator perform well, with the REG estimator performing best. These findings highlight that the choice of estimator may critically depend on the availability and explanatory power of the auxiliary variables.

Table 7 compares the bootstrap variance estimators in terms of %RB and$\%{\mathrm{CP}}_v$. Consistent with the findings for the finite distribution function estimators, the REG estimator shows degraded variance performance at lower quantiles. The IPW estimator maintains coverage close to 95% $\%{\mathrm{CP}}_v$ under Scenario A, but its $\%{\mathrm{CP}}_v$ declined markedly under Scenario B. The DR estimator achieves both low bias and stable $\%{\mathrm{CP}}_v$ across scenarios, indicating reliable variance estimation.

Table 8 compares the quantile estimators in terms of %RB and RRMSE. The REG estimator shows substantial bias at lower quantiles, whereas the IPW estimator performs well under Scenario A but deteriorates under Scenario B. The DR estimator maintains moderate bias and error across both scenarios, yielding comparatively stable performance overall.

Table 9 reports results for the Woodruff confidence intervals of the quantile estimators-$\%{\mathrm{CP}}_{\xi }$, %L, and %U. The IPW estimator attains $\%{\mathrm{CP}}_{\xi }$ close to the nominal 95% under Scenario A, but coverage drops sharply under Scenario B, accompanied by an upward bias in %U, indicating sensitivity to propensity score misspecification. The REG estimator performs well at the middle and upper quantiles, but shows increased %L at lower quantiles. The DR estimator maintains stable $\%{\mathrm{CP}}_{\xi }$ across scenarios, with only a slight upward bias in %U under Scenario B.

6. Conclusions

This study proposed three estimators—Inverse Probability Weighting (IPW), Regression-based estimation (REG), and Doubly Robust estimation (DR)—for reliable estimation of the finite population distribution function and quantiles within a data integration framework that combines probability and nonprobability samples. We examined both theoretical properties and empirical performance. In particular, the DR estimator offers a practical advantage: it retains asymptotic unbiasedness for the finite distribution function provided that either the propensity score model or the outcome regression model is correctly specified, thereby affording robustness to the model misspecification that frequently arises in applied survey settings. Building on this theoretical foundation, we conducted simulation studies using two populations: the synthetic population of [16] and the 2023 Korean Survey of Household Finances and Living Conditions. Across various evaluation metrics, the DR-based procedures showed robust performance, maintaining low relative bias, stable relative root mean squared error, and coverage probabilities close to 95% even when one of the models was misspecified. Notably, DR outperformed IPW and REG when the regression model was inaccurate or the propensity score model was partially misspecified, while also yielding balanced results in the presence of over-representation of high-income households and in lower quantile regions. Furthermore, the composition of auxiliary variables was found to be crucial for estimation performance. Inclusion of covariates with strong explanatory power improved the performance of REG and DR, whereas limited auxiliary information led to increased bias in certain cases. This underscores the importance of selecting appropriate auxiliary variables at the stages of survey design and data integration. Overall, these findings demonstrate that the proposed methods can mitigate the limitations of nonprobability samples and highlight their potential applicability in data environments such as online panel surveys and web-based sources where representativeness is often difficult to achieve.

The main contributions of this study can be summarized in two aspects. First, unlike previous doubly robust (DR) methods that have primarily focused on mean estimation, we extended the approach to the estimation of finite population distribution functions and quantiles. This extension enables more precise and flexible analyses in domains where distributional characteristics such as income, consumption, and health are of central importance. Second, the proposed method enhances the utility of nonprobability samples while being naturally integrated into the framework of probability-based inference, thereby providing an analytical framework well suited for modern survey environments where multiple data sources coexist. Nevertheless, several limitations remain. First, the asymptotic unbiasedness of the DR estimator requires that either the propensity score model or the regression model satisfies certain regularity conditions. When sample sizes are small or the distribution of propensity scores is highly imbalanced, the estimation may become unstable. Second, methodologies for variance estimation in data integration settings are not yet fully established. Conventional bootstrap procedures may overestimate variance, indicating the need for refined theoretical approaches. Third, the present study was conducted under the Missing at Random (MAR) assumption. However, in practice, situations of Not Missing at Random (NMAR) and structural undercoverage occur frequently, highlighting the necessity of developing estimation procedures and diagnostic tools that can address such issues. Future research directions include nonparametric or semiparametric propensity score estimation, integration of high-dimensional auxiliary information through machine learning methods, and applications to a wider range of empirical data sources. Methodological advances along these lines will enable the production of reliable statistics that can accommodate the complexities of real-world survey environments, thereby contributing to evidence-based policymaking using public data.