On valid descriptive inference from non-probability sample

ABSTRACT

We examine the conditions under which descriptive inference can be based directly on the observed distribution in a non-probability sample, under both the super-population and quasi-randomisation modelling approaches. Review of existing estimation methods reveals that the traditional formulation of these conditions may be inadequate due to potential issues of under-coverage or heterogeneous mean beyond the assumed model. We formulate unifying conditions that are applicable to both types of modelling approaches. The difficulties of empirically validating the required conditions are discussed, as well as valid inference approaches using supplementary probability sampling. The key message is that probability sampling may still be necessary in some situations, in order to ensure the validity of descriptive inference, but it can be much less resource-demanding given the presence of a big non-probability sample.

KEYWORDS:

• Non-informative selection
• prediction model
• calibration
• inverse propensity weighting
• sample matching
• model misspecification

1. Introduction

There is a resurgence of interest in the use of non-probability samples. See, for example, Baker et al. (2013) and Elliott and Valliant (2017) for two recent reviews. Such data may arise in situations where probability sampling is either infeasible or too costly. The observations may be obtained from the so-called big-data sources, such as payment transaction data via a specific platform, cellphone call data from a major provider of the service. These big-data non-probability samples can be much larger in size, compared to the more familiar non-probability samples collected from web panel surveys, quota sampling, etc.

Following Rubin (1976) and Little (1982), Smith (1983) considers the so-called super-population (SP) approach to inference from non-probability sample. Under this approach, a prediction model is constructed for the outcome variable of interest, often conditional on some chosen covariates. In particular, Smith (1983) observes an important distinction between analytic and descriptive inference. In analytic inference, the target is the model parameters that are of a theoretical nature; such parameters can never be observed directly no matter how large the sample is. Whereas the targets of descriptive inference are statistics of a given finite population, such that in principle they can be directly observed given a perfect census of the population.

Moreover, Smith (1983) focuses on validity conditions, under which the non-probability sample observation mechanism can be ignored, in the sense that inference can be based on the observed distributions directly, such as the conditional distribution of the outcome variable given the covariates in the sample. The two key validity conditions under the SP approach can be roughly stated as follows: (i) the prediction model is correctly specified for the population units, (ii) the non-probability sample selection mechanism is non-informative, in the sense that the relevant distribution under the population model can be observed in the non-probability sample directly. Similar validity conditions for the SP approach apply in other situations, such as purposive sampling (Royall, 1970), missing data problems (Rubin, 1976).

In this paper, we concentrate on descriptive inference methods that depend on validity conditions in the sense of Smith (1983). Of course, inference is also possible without such validity conditions. For instance, not missing-at-random models (Rubin, 1976) can be used to deal with informative missing data, where the unobserved full-sample outcome distribution is not the same as that among the respondent subsample. Or, the sample likelihood of Pfeffermann, Krieger, and Rinott (1998) can be applied to survey data under informative sampling, where the distribution that holds in the population cannot be directly observed in the sample. See also Pfeffermann (2017) for several other situations where this approach may be relevant. We do not consider such approaches here, which require explicitly modelling the informative observation mechanism of sample selection or measurement.

As reviewed by Elliott and Valliant (2017), there exists another quasi-randomisation (QR) approach to non-probability samples. Under the QR approach, one hypothesises a randomisation model of the non-probability sample inclusion indicator, but treats the outcomes of interest as unknown constants in the population. Though it is clearly inspired by the randomisation approach based on probability sampling, the QR approach is also a model-based approach, based on a model of the sample inclusion indicator instead of a prediction model of the outcome variable under the SP approach. A key motivation is that the correct inclusion probability can be used for any outcome of interest, just like when it is known under probability sampling, whereas the SP approach by nature must be specified differently for different outcome variables. In the context of survey sampling, the QR approach was introduced to deal with nonresponse, where response to survey is modelled as the second phase of selection, in addition to the first phase of sample selection according to a probability sampling design (Oh & Scheuren, 1983).

According to Elliott and Valliant (2017), two key validity conditions are required for the QR approach. (I) The non-probability sample does have a probability sampling mechanism, even though it is unknown. In particular, one assumes that this hypothesised sample inclusion probability is strictly positive for all the population units, so that the only difference to probability sampling is that the inclusion probability is unknown. (II) There exist a set of covariates that ‘fully govern the sampling mechanism’. In other words, the sample inclusion probability is a function of these covariates.

Thus there are two model-based approaches to inference from non-probability sample. Under the SP approach, one models the outcome variable conditional on the realised sample inclusion indicators, whereas under the QR approach, one models the sample inclusion indicators, but treats the outcomes as unknown constants. Although one may envisage the outcomes as the realised values of random variables, a fully specified model of the outcome variable will not be required under the QR approach, given suitable validity conditions. Similarly, although one acknowledges that the sample selection mechanism may be critical to the SP approach, a fully specified model of the inclusion indicator will not be required under the SP approach, given suitable validity conditions.

It is possible to construct estimators that combine both the models of outcome and sample inclusion indicator, in a manner such that the estimator is consistent as long as one of the two models hold. Over the recent years, it is becoming common to refer to this estimation approach as ‘doubly robust’ (Kang & Schafer, 2007; Kim & Haziza, 2014; Robins, Rotnitzky, & Zhao, 1994). Notice that the traditional generalised regression estimator in survey sampling is doubly robust in the same sense, except that here the randomisation mechanism is actually known. Nevertheless, it is a fact that in the debate between model-based and design-based inference from probability sampling, either side questioned the ‘robustness’ of the other.

The rest of the paper is organised as follows. In Section 2, we review the estimation methods from non-probability sample which do require validity conditions. Although these have been roughly stated above, a closer examination under both modelling perspectives reveals nuances across the different estimators. Moreover, we shall highlight the potential challenges of under-coverage and heterogeneous means beyond the assumed model. The traditional formulation of validity conditions is inadequate in both regards. We outline a set of unified validity conditions in Section 3, which are formulated non-parametrically and encompasses both the modelling approaches. Post-stratification and calibration estimators are considered in light of these conditions. However, as will be discussed, a key difficulty in practice is that the validity conditions may be impossible to verify empirically based only on the data used for the estimation. Finally, we outline shortly in Section 4 some valid approaches given a supplementary probability sampling of the outcome of interest, followed by a brief summary in Section 5.

The key message is that probability sampling may still be necessary in some situations, in order to ensure the validity of descriptive inference, but it can be much less resource-demanding given the presence of a big non-probability sample. In fact, the bigger the non-probability sample, the better it is.

2. Review of existing approaches

Denote by U the population of known size N. Let each population unit be associated with an outcome of interest, denoted by $y_i$, for $i\in U$. Denote by B the observed nonprobability sample of size $n_B$. A common assumption to all the estimators we discuss below is that B does not contain any out-of-scope units, such that $B\subset U$, and there are no duplicated units in B. Let ${\delta }_i=1$ if $i\in B$, and 0 if $i\in U\setminus B$. Let $y_i$ be observed for all the units in B, and let $y_B=\{y_i;i\in B\}$. To fix the idea, let $Y=\sum _{i\in U}{y}_i$ be the population total that is the target of descriptive inference. Let $x_B=\{x_i;i\in B\}$ in cases where any relevant covariates $x_i$ are available in the sample B. Let $X=\sum _{i\in U}{x}_i$ be the population totals and let $\overline{X}=X/N$. Given $x_B$, one can have two situations depending on whether $(X,\overline{X})$ are known or not. In the case they are unknown, it may still be possible that there exists a second probability sample S, for $S\subset U$, in which $x_i$ is observed, so that $(X,\overline{X})$ can be estimated based on the sample S.

2.1. B-sample expansion estimator

Consider first the most basic situation where only $y_B$ is observed, and no relevant covariates are available at all. Let ${\overline{y}}_B=\sum _{i\in B}{y}_i/n_B$ be the B-sample mean. The B-sample expansion estimator of Y is given by (1) $\hat{Y}=N{\overline{y}}_B.$(1) Under the SP approach, let ${\mu }_i=E(y_i|{\delta }_i,i\in U)$ be the conditional expectation of $y_i$ given ${\delta }_i$, for any $i\in U$, where both ${\delta }_i$ and $y_i$ are treated as random variables. Provided the conditional expectation is the same as the unconditional expectation given either ${\delta }_i=1$ or ${\delta }_i=0$, for any $i\in U$, denoted by (2) $\mu =\mu ({\delta }_i=1)=E(y_i|i\in B)=E(y_i;i\in U),$(2) we have $E({\overline{y}}_B-Y/N|B)=\sum _{i\in B}\mu /n_B-\mu =0$ such that the B-sample expansion estimator is prediction unbiased for Y. We shall refer to (2(2) $\mu =\mu ({\delta }_i=1)=E(y_i|i\in B)=E(y_i;i\in U),$(2) ) as the SP assumption, which is a validity condition for the B-sample expansion estimator under the SP approach.

Under the QR approach, where $y_i$ is treated as a fixed constant, let $p_i=\text{Pr}({\delta }_i=1;y_i,i\in U)$ be the inclusion probability of any population unit that is associated with the value $y_i$. The notation ‘;’ is used here instead of ‘$|$’ because, strictly speaking, $p_i$ is not a conditional probability now that $y_i$ is not conceived as the realised value of a random variable under the QR approach. Now, provided the inclusion probability is the same for any $i\in U$, (3) $p_i=p$(3) we have $\tilde{Y}=\sum _{i\in B}{y}_i/p$ is unbiased for Y, since $\begin{array}{rl}E\left(\sum _{i\in B}{y}_i/p\right)& =\sum _{i\in U}E({\delta }_i;y_i,i\in U)y_i/p\\ & =\sum _{i\in U}p{y}_i/p=Y.\end{array}$ In reality, p is unknown. It is natural to estimate it by $\hat{p}=n_B/N$ under (3(3) $p_i=p$(3) ), which yields (1(1) $\hat{Y}=N{\overline{y}}_B.$(1) ) as the resulting plug-in estimator. It follows that the QR assumption (3(3) $p_i=p$(3) ) is the key validity condition, which ensures that the B-sample expansion estimator is consistent for Y, as $N\to \mathrm{\infty }$ and $n_B/N=O_p\left(1\right)$ asympotically.

In summary, the B-sample expansion estimator (1(1) $\hat{Y}=N{\overline{y}}_B.$(1) ) can be motivated under both the SP and QR approaches, given validity conditions (2(2) $\mu =\mu ({\delta }_i=1)=E(y_i|i\in B)=E(y_i;i\in U),$(2) ) and (3(3) $p_i=p$(3) ), respectively.

2.2. B-sample calibration estimator

Suppose relevant covariates $x_B$ are available in the sample B. The population totals X may be either known or unknown. In the latter case, suppose they can be estimated from a second probability sample S. The B-sample calibration estimator of Y is given by (4) $\begin{array}{rl}\hat{Y}& =\sum _{i\in B}{w}_i{y}_i\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\\ & \left\{\begin{array}{ll}\sum _{i\in B}{w}_i{x}_i=X& \mathrm{i}\mathrm{f}\mathrm{k}\mathrm{n}\mathrm{o}\mathrm{w}\mathrm{n}X,\\ \sum _{i\in B}{w}_i{x}_i=\hat{X}\left(S\right)& \mathrm{i}\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{k}\mathrm{n}\mathrm{o}\mathrm{w}\mathrm{n}X,\end{array}\right.\end{array}$(4) where $\hat{X}\left(S\right)$ is some consistent S-sample estimator, as the S-sample size increases, and the weights $w_B=\{w_i;i\in B\}$ are calibrated in a way depending on the availability of X.

To actually compute the estimator (4(4) $\begin{array}{rl}\hat{Y}& =\sum _{i\in B}{w}_i{y}_i\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\\ & \left\{\begin{array}{ll}\sum _{i\in B}{w}_i{x}_i=X& \mathrm{i}\mathrm{f}\mathrm{k}\mathrm{n}\mathrm{o}\mathrm{w}\mathrm{n}X,\\ \sum _{i\in B}{w}_i{x}_i=\hat{X}\left(S\right)& \mathrm{i}\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{k}\mathrm{n}\mathrm{o}\mathrm{w}\mathrm{n}X,\end{array}\right.\end{array}$(4) ), one needs to choose a set of initial weights, denoted by $a_B=\{a_i;i\in B\}$, and a distance function such as $\sum _{i\in B}(w_i-a_i{)}^2/a_i$ between the initial and calibrated weights (Deville & Särndal, 1992). In the case of (5) $a_i=1/p_i,$(5) where $p_i$ is the true B-sample inclusion probability, for $p_i>0$, the calibration estimator is consistent, as $N\to \mathrm{\infty }$ and $n_B/N=O_p\left(1\right)$, given mild regularity conditions in addition. However, insofar as one cannot manage to set the initial weights (5(5) $a_i=1/p_i,$(5) ), the calibration estimator is unmotivated from the QR perspective.

Next, under the SP approach, suppose the SP${}_x$ assumption given by (6) $E(y_i|x_i,i\in U)=\mu \left(x_i\right)=x_i^{\mathrm{\top }}\text{β},$(6) which relates the conditional expectation of $y_i$ linearly to the given $x_i$, and (7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) by which the B-sample selection is non-informative given $x_i$. We have then $E\left(\sum _{i\in B}{w}_i{y}_i-Y|x_U\right)=E\left(\sum _{i\in B}{w}_i{x}_i^{\mathrm{\top }}\text{β}\right)-X^{\mathrm{\top }}\text{β}=0$ given $\sum _{i\in B}{w}_i{x}_i=X$, regardless of the initial weights $a_B$. Otherwise, this expectation would tend to 0, provided $\hat{X}\left(S\right)$ is an asymptotically unbiased estimator of X, under some suitable asymptotic setting. It follows that the assumptions (6(6) $E(y_i|x_i,i\in U)=\mu \left(x_i\right)=x_i^{\mathrm{\top }}\text{β},$(6) ) and (7(7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) ) are the key validity conditions for the B-sample calibration estimator under the SP approach.

The estimator (4(4) $\begin{array}{rl}\hat{Y}& =\sum _{i\in B}{w}_i{y}_i\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\\ & \left\{\begin{array}{ll}\sum _{i\in B}{w}_i{x}_i=X& \mathrm{i}\mathrm{f}\mathrm{k}\mathrm{n}\mathrm{o}\mathrm{w}\mathrm{n}X,\\ \sum _{i\in B}{w}_i{x}_i=\hat{X}\left(S\right)& \mathrm{i}\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{k}\mathrm{n}\mathrm{o}\mathrm{w}\mathrm{n}X,\end{array}\right.\end{array}$(4) ) becomes the B-sample post-stratification estimator in the special case where $x_i$ is the post-stratum dummy index. For the QR approach, one can set $a_i$ to be the inverse post-stratum B-sample fraction, which is equivalent to introducing the QR assumption (3(3) $p_i=p$(3) ) in each post-stratum separately. This QR${}_x$ assumption provides then a validity condition for the B-sample post-stratification estimator under the QR approach. For the SP approach, the two assumptions (6(6) $E(y_i|x_i,i\in U)=\mu \left(x_i\right)=x_i^{\mathrm{\top }}\text{β},$(6) ) and (7(7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) ) remain formally the same.

2.3. B-sample inverse propensity weighting

Suppose relevant covariates $x_B$ are available in the sample B. The B-sample inverse propensity weighting (IPW) estimator is constructed under the QR approach. Suppose (8) $p_i=p(x_i;\eta )>0$(8) i.e., the B-sample inclusion probability is completely determined given $x_i$, in the strictly positive parametric form $p(x_i;\eta )$, which may as well be referred to as the QR${}_x$ assumption. Provided $x_i$ is known for all the units in the population, η can be estimated, say, by a population estimating equation $\sum _{i\in U}H({\delta }_i;\eta )=0,$ where $E\left[H\right({\delta }_i;\eta \left)\right]=0$. Otherwise, suppose $x_S$ is observed in a second probability sample S, one can use the pseudo population estimating equation $\sum _{i\in S}{d}_{i}H({\delta }_i;\eta )=0$ (Kim & Wang, 2018), where $d_i$ is the sampling weight, for $i\in S$, or some S-sampling design-consistent adjustment of it. This requires that one is able to observe ${\delta }_i$ for each unit i in S, in other words the two samples S and B can be matched, which is an important assumption in terms of application. To ensure that $H({\delta }_i;\eta )$ is the same in both of these two estimating equations, i.e., whether $i\in S$ or just $i\in U$, one needs to assume that S-sampling from U is non-informative for ${\delta }_i$, so that (9) $\text{Pr}({\delta }_i=1|x_i,i\in S)=\text{Pr}({\delta }_i=1|x_i,i\in U).$(9) Notice that, given non-informativeness (9(9) $\text{Pr}({\delta }_i=1|x_i,i\in S)=\text{Pr}({\delta }_i=1|x_i,i\in U).$(9) ), we have $E\left[H\right({\delta }_i;\eta \left)\right]=0$ for all $i\in s$, such that one can also use the unweighted S-sample estimating equation, which is given by $\sum _{i\in S}H({\delta }_i;\eta )=0$ instead of the pseudo population estimating equation. Having obtained the parameter estimate $\hat{\eta }$, one obtains ${\hat{p}}_i=p(x_i;\hat{\eta })$ and the B-sample IPW estimator (10) $\hat{Y}=\sum _{i\in B}{y}_i/{\hat{p}}_i,$(10) which is consistent for Y under mild regularity conditions, if $\hat{\eta }$ is consistent for η under some suitable asymptotic setting. It follows that the QR${}_x$ assumption (8(8) $p_i=p(x_i;\eta )>0$(8) ) is its key validity condition, whereas the non-informativeness assumption (9(9) $\text{Pr}({\delta }_i=1|x_i,i\in S)=\text{Pr}({\delta }_i=1|x_i,i\in U).$(9) ) is needed in addition when $x_i$ is only available in a probability sample S instead of the population.

2.4. Another B-sample IPW estimator

Elliott and Valliant (2017) discuss another IPW estimator (10(10) $\hat{Y}=\sum _{i\in B}{y}_i/{\hat{p}}_i,$(10) ), where $p_i$ is obtained with the help of a second so-called reference probability sample S, and is given by (11) $p_i\propto \text{Pr}(S_i=1|x_i,i\in U)\frac{\text{Pr}({\delta }_i=1|x_i,i\in B\cup S)}{\text{Pr}(S_i=1|x_i,i\in B\cup S)},$(11) where $S_i=1$ if $i\in S$ and 0 if $i\in U\setminus S$, and to fix the idea one may suppose $S\cap B=\mathrm{\varnothing }$. First, the QR${}_x$ assumption (8(8) $p_i=p(x_i;\eta )>0$(8) ) is retained. The definition of $p_i$ by (11(11) $p_i\propto \text{Pr}(S_i=1|x_i,i\in U)\frac{\text{Pr}({\delta }_i=1|x_i,i\in B\cup S)}{\text{Pr}(S_i=1|x_i,i\in B\cup S)},$(11) ) can then be motivated as follows: $\begin{array}{rl}& \frac{\text{Pr}({\delta }_i=1|x_i,i\in U)}{\text{Pr}(S_i=1|x_i,i\in U)}\\ & \propto \frac{\text{Pr}(x_i|{\delta }_i=1,i\in U)}{\text{Pr}(x_i|S_i=1,i\in U)}\\ & \left[\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}.\mathrm{t}\mathrm{o}\frac{\text{Pr}({\delta }_i=1|i\in U)}{\text{Pr}(S_i=1|i\in U)}\right]\\ & \propto \frac{\text{Pr}(x_i|{\delta }_i=1,i\in B\cup S)}{\text{Pr}(x_i|S_i=1,i\in B\cup S)}\\ & \propto \frac{\text{Pr}({\delta }_i=1|x_i,i\in B\cup S)}{\text{Pr}(S_i=1|x_i,i\in B\cup S)}\\ & \left[\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}.\mathrm{t}\mathrm{o}\frac{\text{Pr}({\delta }_i=1|i\in B\cup S)}{\text{Pr}(S_i=1|i\in B\cup S)}\right]\end{array}$ provided the S-sample inclusion probability is also fully determined by $x_i$ in the sense of (8(8) $p_i=p(x_i;\eta )>0$(8) ). Thus the validity condition for the IPW estimator (10(10) $\hat{Y}=\sum _{i\in B}{y}_i/{\hat{p}}_i,$(10) ) based on (11(11) $p_i\propto \text{Pr}(S_i=1|x_i,i\in U)\frac{\text{Pr}({\delta }_i=1|x_i,i\in B\cup S)}{\text{Pr}(S_i=1|x_i,i\in B\cup S)},$(11) ) is that the QR${}_x$ assumption (8(8) $p_i=p(x_i;\eta )>0$(8) ) holds for both the samples, given the same $x_i$.

We make two observations. First, despite the superficial resemblance to the propensity scoring method of Rosenbaum and Rubin (1983), the above argument for $p_i$ is not the same. As Rosenbaum and Rubin (1983) state clearly before their first enumerated equation, ‘in this paper, the N units in the study are viewed as a simple random sample from some population’, where N is the size of the combined sample of treatment and non-treatment. The analogy to this combined sample is $B\cup S$ here. However, it is generally untenable that $B\cup S$ can be treated as a simple random sample from the population. Second, for any given probability sample S, it is possible to identify the variables that determine the designed inclusion probability, denoted by ${\pi }_i=\pi \left(z_i\right)$, for $i\in U$. There arises thus a question, ‘what if $\pi \left(z_i\right)$ differs considerably from $p(x_i,\hat{\eta })$?’ Moreover, one may have more than one probability sample in which $x_i$ is observed. There arises then a question, ‘which reference sample should one use?’

2.5. Sample matching estimator

Rivers (2007) applies the SP approach in situations where a second probability sample S is available. Yang and Kim (2018) study mass imputation methods, which include the matching estimator of Rivers (2007) as a special case. The sample matching (SM) estimator is given by (12) $\hat{Y}=\sum _{i\in S}{d}_i{\hat{y}}_i,$(12) where ${\hat{y}}_i=y_{k_i}$, for $k_i=\arg \underset{j\in B}{min}\parallel x_i-x_j\parallel$ based on a chosen metric $\parallel \cdot \parallel$. That is, $y_{k_i}$ is the nearest-neighbour (NN) imputation value from the B-sample for $i\in S$.

To focus on the idea, assume for the moment exact matching is the case, where $x_{k_i}=x_i=x$ for all $i\in S$ and $k_i\in B$. We have then $E({\hat{y}}_i|x_i=x)=E(y_{k_i}|x_{k_i}=x,k_i\in B)$, which is the same as $E(y_i|x_i=x,i\in B)$ as if the unit i were in B. Given the non-informativeness assumption (7(7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) ) for the B-sample, which Yang and Kim (2018) call the ‘ignorability’ assumption, we have $\begin{array}{rl}& E\left[\sum _{i\in S}{d}_{i}E\left({\hat{y}}_i|x_i\right)\right]\\ & =E\left[\sum _{i\in S}{d}_{i}E(y_i|x_i,i\in B)\right]\\ & =E\left[\sum _{i\in S}{d}_{i}E(y_i|x_i,i\in U)\right]\\ & =\sum _{i\in U}E(y_i|x_i,i\in U)=E\left(Y|x_U\right).\end{array}$ With respect to both the population model and the design of S, the SM estimator (12(12) $\hat{Y}=\sum _{i\in S}{d}_i{\hat{y}}_i,$(12) ) is prediction unbiased for Y. Notice that in the case of S = U, the SM estimator is just an NN-imputation method. Whether S = U or not, the NN-imputed SM estimator is likely to be less efficient than a prediction-imputed SM estimator $\hat{Y}=\sum _{i\in S}{d}_{i}E(y_i|x_i;\widehat{\text{β}}\left(B\right))$ whenever a correct parametric specification of the conditional mean (via β) is possible. The simulation results of Yang and Kim (2018) show that NN-imputation is less efficient than imputation based on semi-parametric generalised additive models.

Now, it is not difficult to see that the consistency of the SM estimator (12(12) $\hat{Y}=\sum _{i\in S}{d}_i{\hat{y}}_i,$(12) ) can be established, given asymptotic exact matching instead, i.e., (13) $\parallel x_i-x_{k_i}\parallel \to 0\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y},$(13) for any $i\in S$, as $N\to \mathrm{\infty }$ and $n_B/N=O_p\left(1\right)$. Yang and Kim (2018) make the assumption of ‘common support’ to the same effect. To ensure that $E(y_i|x_i,i\in U)$ does not change abruptly as the value $x_i$ varies, Yang and Kim (2018) assume that $E(y_i|x_i,i\in U)$ is continuous differentiable. Or, one may adopt the SP${}_x$ assumption below: (14) $\begin{array}{rl}& \parallel E(y_i|x_i,i\in U)-E(y_j|x_j,j\in U)\parallel \\ & =O(\parallel x_i-x_j\parallel )\mathrm{a}\mathrm{s}N\to \mathrm{\infty }\end{array}$(14) (Chen & Shao, 2000, Theorem 1). It follows that the assumptions (7(7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) ), (13(13) $\parallel x_i-x_{k_i}\parallel \to 0\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y},$(13) ) and (14(14) $\begin{array}{rl}& \parallel E(y_i|x_i,i\in U)-E(y_j|x_j,j\in U)\parallel \\ & =O(\parallel x_i-x_j\parallel )\mathrm{a}\mathrm{s}N\to \mathrm{\infty }\end{array}$(14) ) are the key validity conditions for the consistency of the SM estimator (12(12) $\hat{Y}=\sum _{i\in S}{d}_i{\hat{y}}_i,$(12) ).

We make two observations. First, an attractive feature of the NN-imputation is that the imputed sample S looks more realistic and natural than, say, by the regression prediction imputation. However, unless the S-sampling is non-informative, the NN-imputed S-sample will not resemble the true S-sample that could have been observed, since $E({\hat{y}}_i|x_i,i\in S)=E(y_i|x_i,i\in U)\ne E(y_i|x_i,i\in S),$ where the inequality is the case unless S-sampling is non-informative in the sense of (7(7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) ). Second, for any other covariate $z_i\ne x_i$, including when $z_i$ contains the S-sample design variables, we have $\begin{array}{rl}& E({\hat{y}}_i|z_i,x_i,i\in S)=E(y_i|x_i,i\in U)\\ & \ne E(y_i|z_i,x_i,i\in U)\end{array}$ unless $y_i$ and $z_i$ are conditionally independent of each other given $x_i$. This is a general problem for statistical matching of variables associated with distinct units, i.e., $y_i$ associated with $x_i$ for some $i\in B$ and $z_i$ associated with the same value $x_i$ but for some different unit in S. The following example illustrates both remarks above.

Example

Let $y_i$ be independent of $x_i\sim \text{Unif}(0,1)$, for any $i\in U$. Then, the SP${}_x$ assumption (14(14) $\begin{array}{rl}& \parallel E(y_i|x_i,i\in U)-E(y_j|x_j,j\in U)\parallel \\ & =O(\parallel x_i-x_j\parallel )\mathrm{a}\mathrm{s}N\to \mathrm{\infty }\end{array}$(14) ) holds trivially, as long as the marginal expectation $E\left(y_i\right)$ exists. Next, suppose simple random sample B, so that the non-informative assumption (7(7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) ) holds, and $E({\hat{y}}_i|x_i,i\in S)=E(y_i|i\in U)$ regardless of the exact matching assumption. Suppose stratified simple random S-sampling with two strata of different sampling fractions, so that the S-sample inclusion probability is not a constant. Then, the S-sampling is informative (given $x_i$) as long as the population stratum means are different, since $E({\overline{y}}_S|x_S,S)=E\left({\overline{y}}_S|S\right)\ne E\left(\overline{Y}\right)=E\left(\overline{Y}|x_U\right),$ where ${\overline{y}}_S$ is the true S-sample mean that is unknown, since $y_i$ is not observed in S. It follows that the NN-imputed S-sample $\{{\hat{y}}_i;i\in S\}$ would look like a sample generated by simple random sampling, rather than the actual stratified sampling. Moreover, the SM-estimator of stratum means, corresponding to say $z_i=1,2$, respectively, will be biased for the population stratum means.

3. More generally on validity conditions

Non-informative selection in form of (7(7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) ) or (9(9) $\text{Pr}({\delta }_i=1|x_i,i\in S)=\text{Pr}({\delta }_i=1|x_i,i\in U).$(9) ) is a critical condition for all the methods in Section 2, which make use of auxiliary variable $x_i$. Two kinds of possible violation of these assumptions are worth noting.

First, Kim and Rao (2018) point out an important issue that has not received sufficient attention in these methods, namely B-sample under-coverage is the case if some population units have in fact zero chance of being included in it. Under the SP approach, extrapolation of the conditional distribution of $y_i$ in the B-sample to these population units can only be based on subjective beliefs but not empirical evidence. The QR approach is equally affected, since randomisation inference would have been invalidated even if $p_i$ were known for all the B-sample units, let alone when it is unknown and needs to be estimated. To address the issue, Kim and Rao (2018) consider a two-phase SM estimator. Let the S-sample be partitioned into $S_1$ and $S_0$, such that $S_1=\{i;p_i>0\}$ and $S_0=\{i;p_i=0\}$. First, estimate this unobserved partition via the B-sample support: ${\hat{S}}_1=\left\{i;\underset{j\in B}{min}\parallel x_i-x_j\parallel <ϵ\right\}.$ Each S-sample unit that is unsupported in the B-sample ε-neighbourhood is assigned to ${\hat{S}}_0$. Let us suppose this partition estimator is consistent in the following sense: $|{\hat{S}}_1\cup S_1|/|{\hat{S}}_1\cap S_1|\to 1\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y},$ as $N\to \mathrm{\infty }$ and $ϵ\to 0$. Next, the two-phase SM estimator is given as $\hat{Y}=\sum _{i\in {\hat{S}}_1}{d}_i{w}_{2i}{\hat{y}}_i,$ where $\sum _{i\in {\hat{S}}_1}{d}_i{w}_{2i}{x}_i=\sum _{i\in S}{d}_i{x}_i$. In other words, the under-coverage is dealt with by the calibration of the weights $w_{2i}$. This can be motivated, provided the conditional mean $E(y_i|x_i,p_i=0)$ can be linearly related to $x_i$, and the relationship is the same for the units with $p_i>0$, i.e., the under-coverage is non-informative for the SP linear model.

Second, insofar one requires either an assumption of SP${}_x$ (7(7) $E(y_i|x_i,i\in U)=E(y_i|x_i,i\in B)$(7) ) or QR${}_x$ (9(9) $\text{Pr}({\delta }_i=1|x_i,i\in S)=\text{Pr}({\delta }_i=1|x_i,i\in U).$(9) ), there is always the possibility of heterogeneous mean, beyond what is controlled by the chosen $x_i$. Let $U_x=\{i;x_i=x,i\in U\}$ be of the size $N_x$. Under the SP approach, which models the mean ${\mu }_i$ of unit i by $\mu \left(x_i\right)$, heterogeneous mean is the case if ${\mu }_i\ne \mu \left(x_i\right)$, despite (15) $\mu \left(x\right)=\sum _{i\in U_x}{\mu }_i/N_x$(15) and $\mu \left(x_i\right)$ is statistically correct in that the ${\mu }_i$'s average to $\mu \left(x\right)$ over all the units in $U_x$. Under the QR approach, heterogeneous mean is the case if $p_i\ne p\left(x_i\right)$, despite (16) $p\left(x\right)=\sum _{i\in U_x}{p}_i/N_x.$(16) Let us illustrate the concept of heterogeneous mean with a simple example.

Example

Let $x\equiv 1$, such that $\mu \left(x_i\right)=\mu$, for all $i\in U$. Let $U=U_1\cup U_0$ be a partition. Let $U_1$ be of size $N_1$ and with mean ${\mu }_i=\mu \left(1\right)$, for all $i\in U_1$; let $U_0$ be of size $N_0$ and with mean ${\mu }_i=\mu \left(0\right)$, for all $i\in U_0$. Suppose $\mu \left(1\right)\ne \mu \left(0\right)$. Let $\mu =\mu \left(1\right)N_1/N+\mu \left(0\right)N_0/N$. Then, ${\mu }_i\ne \mu$ for any $i\in U$, but we still have $\sum _{i\in U}{\mu }_i/N=\mu$, satisfying (15(15) $\mu \left(x\right)=\sum _{i\in U_x}{\mu }_i/N_x$(15) ).

Heterogeneous mean affects the SP and QR approaches differently. Given (15(15) $\mu \left(x\right)=\sum _{i\in U_x}{\mu }_i/N_x$(15) ), assuming ${\mu }_i=\mu \left(x\right)$ for $i\in U_x$ is prediction unbiased, despite heterogeneous mean, since $\sum _{i\in U_x}\left[E\right(y_i|{\delta }_i)-\mu (x\left)\right]=\sum _{i\in U_x}[{\mu }_i-\mu (x\left)\right]=0.$ Meanwhile, given (16(16) $p\left(x\right)=\sum _{i\in U_x}{p}_i/N_x.$(16) ), assuming $p_i=p\left(x\right)$ for $i\in U_x$ yields $\begin{array}{rl}& E\left(\sum _{i\in U_x}\frac{{\delta }_i{y}_i}{p\left(x\right)}\right)-\sum _{i\in U_x}{y}_i\\ & =p(x{)}^{-1}\sum _{i\in U_x}(p_i-p(x))y_i\ne 0\end{array}$ in which case the IPW estimator under the QP approach may be biased, despite the model of $p_i$ is statistically correct in the sense of (16(16) $p\left(x\right)=\sum _{i\in U_x}{p}_i/N_x.$(16) ).

The discussion above suggests that the formulation of validity conditions in Section 2 is inadequate in the presence of under-coverage and mean heterogeneity. Below we first reformulate the validity conditions, which cover both the SP and QR approaches, despite the presence of under-coverage and mean heterogeneity. We elaborate and illustrate these conditions for the post-stratification and calibration estimators. Finally, we discuss the difficulties of verifying these validity conditions empirically.

3.1. Non-parametric asymptotic (NPA) non-informativeness

We start by noticing that the B-sample mean equals to the population mean, denoted by ${\overline{y}}_B=\overline{Y}$, provided $\begin{array}{rl}& {\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)\\ & =\frac{1}{N}\sum _{i\in U}{\delta }_i{y}_i-\left(\frac{1}{N}\sum _{i\in U}{\delta }_i\right)\left(\frac{1}{N}\sum _{i\in U}{y}_i\right)=0,\\ & E_N\left({\delta }_i\right)=\sum _{i\in U}{\delta }_i/N>0,\end{array}$ where $E_N$ and ${\mathrm{C}\mathrm{o}\mathrm{v}}_N$ denote, respectively, expectation and covariance with respect to the empirical distribution function that places point mass 1/N on each population unit. This provides an empirical formulation of the non-informativeness of the B-sample observation mechanism with respect to the outcome of interest. Similar expressions have appeared in various discussions of the potential sample mean bias due to the observation mechanism, such as unequal probability sampling (Rao, 1966), survey nonresponse (Bethlehem, 1988) or big data (Meng, 2018). It motivates the following non-parametric asymptotic (NPA) non-informativeness assumption in the absence of any covariates:

(17) $\begin{array}{rl}\underset{N\to \mathrm{\infty }}{lim}{\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)=0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ \underset{N\to \mathrm{\infty }}{lim}{E}_N\left({\delta }_i\right)=p>0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}$(17) The NPA assumption (17(17) $\begin{array}{rl}\underset{N\to \mathrm{\infty }}{lim}{\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)=0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ \underset{N\to \mathrm{\infty }}{lim}{E}_N\left({\delta }_i\right)=p>0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}$(17) ) encompasses both the SP and QR approaches. For the SP approach, taking the conditional expectation of $y_i$'s conditional on the ${\delta }_i$'s yields $\begin{array}{rl}& E({\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)|{\delta }_U)\\ & =\frac{1}{N}\sum _{i\in U}{\delta }_i{\mu }_i-\left(\frac{1}{N}\sum _{i\in U}{\delta }_i\right)\left(\frac{1}{N}\sum _{i\in U}{\mu }_i\right)\to 0\end{array}$ given NPA non-informative B-selection, where $\sum _{i\in U}{\delta }_i/N>0$ given non-negligible B-selection in addition. Under this condition, the B-sample expansion estimator (1(1) $\hat{Y}=N{\overline{y}}_B.$(1) ) is asymptotically prediction unbiased from the SP perspective. For the QR approach, taking the expectation of ${\delta }_i$ s with the $y_i$ s being constants yields $\begin{array}{rl}& E({\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i);y_U)\\ & =\frac{1}{N}\sum _{i\in U}{p}_i{y}_i-\left(\frac{1}{N}\sum _{i\in U}{p}_i\right)\left(\frac{1}{N}\sum _{i\in U}{y}_i\right)\to 0,\\ & E(E_N\left({\delta }_i\right))=\sum _{i\in U}{p}_i/N\to p>0.\end{array}$ In particular, the NPA assumption (17(17) $\begin{array}{rl}\underset{N\to \mathrm{\infty }}{lim}{\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)=0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ \underset{N\to \mathrm{\infty }}{lim}{E}_N\left({\delta }_i\right)=p>0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}$(17) ) allows for $0\le p_i\le 1$, so that the B-sample expansion estimator (1(1) $\hat{Y}=N{\overline{y}}_B.$(1) ) remains consistent from the QR perspective, even in the presence of under-coverage of the units with $p_i=0$ or non-representative units with $p_i=1$.

Example

Let $U=U_1\cup U_0$ be a partition. Let $U_1$ be of size $N_1$, where $p_i\equiv 1$ for $i\in U_1$; let $U_0$ be of size $N_0$, where $p_i\equiv 0$ for $i\in U_0$. Despite under-coverage of $B\equiv U_1$, the first NPA condition implies ${\overline{y}}_B-\overline{Y}\to 0$, given the second condition $N_1/N\to p>0$.

3.2. Post-stratification estimator

Consider post-stratification by $x_i$, for $i\in U$. Provided the assumption (17(17) $\begin{array}{rl}\underset{N\to \mathrm{\infty }}{lim}{\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)=0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ \underset{N\to \mathrm{\infty }}{lim}{E}_N\left({\delta }_i\right)=p>0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}$(17) ) holds within each post-stratum, the B-sample post-stratification estimator is asymptotically unbiased from both the SP and QR perspectives. Below we consider the QR approach. The SP approach is a special case of the calibration estimator discussed in Section 3.3.

Consider first the hypothetical estimator with known $p_x=\sum _{i\in U_x}{p}_i/N_x$: $\tilde{Y}=\sum _x\sum _{i\in U_x}{\delta }_i{y}_i/p_x.$ To fix the idea for variance estimation, suppose independent Bernoulli distribution of ${\delta }_i$ with probability $p_i$, where $0\le p_i\le 1$. The variance of $\tilde{Y}$ is then given by $V\left(\tilde{Y}\right)=\sum _x\sum _{i\in U_x}{p}_i{y}_i^2/p_x^2-\sum _x\sum _{i\in U_x}{p}_i^2{y}_i^2/p_x^2.$ An unbiased estimator of the first term of the variance, denoted by ${\tau }_1$ is given by ${\hat{\tau }}_1=\sum _x\sum _{i\in U_x}{\delta }_i{y}_i^2/p_x^2=\sum _x{p}_x^{-2}\sum _{i\in B_x}{y}_i^2,$ where $B_x=B\cap U_x$. An unbiased estimator of the second term, denoted by ${\tau }_2$ is given by $\begin{array}{rl}{\hat{\tau }}_2& =\sum _x{p}_x^{-2}\sum _{i\in U_x}{\delta }_i{p}_i{y}_i^2\\ & =\sum _x{p}_x^{-1}\sum _{i\in U_x}{\delta }_i{y}_i^2=\sum _x{p}_x^{-1}\sum _{i\in B_x}{y}_i^2,\end{array}$ where the second equality follows given the additional QR${}_x$ assumption, i.e., $p_i=p_x$ for $i\in U_x$. Putting ${\hat{\tau }}_1$ and ${\hat{\tau }}_2$ together, we obtain $\hat{V}\left(\tilde{Y}\right)=\sum _x(p_x^{-1}-1)p_x^{-1}\sum _{i\in B_x}{y}_i^2.$ Now, the post-stratification estimator, denoted by $\hat{Y}$, is obtained from $\tilde{Y}$ on replacing $p_x$ by ${\hat{p}}_x=n_{xB}/N_x$, where $n_{xB}$ is the observed size of $B_x$ and $N_x$ is the known post-stratum population size. Expanding ${\hat{p}}_x$ around $p_x$ (i.e., linearisation) would yield an asymptotically valid estimator of the unconditional variance of $\hat{Y}$.

3.3. Calibration estimator

The post-stratification estimator is infeasible, in cases when the B-sample contains empty cells, or when the population size $N_x$ is not all known. Let $\begin{array}{rl}{t}_i& =(t_{1i},t_{2i},\dots t_{Ki}{)}^{\mathrm{\top }}\\ & =(t_1\left(x_i\right),t_2\left(x_i\right),\dots t_K\left(x_i\right){)}^{\mathrm{\top }}=t\left(x_i\right)\end{array}$ be a vector of many-to-one mappings of $x_i$, such that the population total $T=\sum _{i\in U}{t}_i$ is known, and the sample total $t=\sum _{i\in B}{t}_i$ has only non-zero components.

As discussed for the calibration estimator in Section 2, generally one is not able to set the initial weight to be the inverse of B-sample inclusion probability in practice. Suppose one simply starts with initial equal weights $a_i=N/n_B$ for all $i\in B$. The linear calibration estimator (Deville & Särndal, 1992) is given by $\hat{Y}=\sum _{i\in B}{w}_i{y}_i,$ where the weights $\{w_i;i\in B\}$ minimise the distance to $\{a_i;i\in B\}$ as measured by $\begin{array}{rl}& \sum _{i\in B}(w_i-N/n_B{)}^2\\ & =\sum _t\left(\sum _{i\in B_t}{w}_i^2-2\left(N/n_B\right)\sum _{i\in B_t}{w}_i+n_{tB}(N/n_{tB}{)}^2\right)\end{array}$ subjected to the constraints $\sum _{i\in B}{w}_i{t}_i=T$, where $B_t=\{i;t_i=t,i\in B\}$ and $n_{tB}>0$. It follows that $w_i=w_t$, for $i\in B_t$, since the only thing that matters to the calibration constraints is $\sum _{i\in B_t}{w}_i$ now that $t_i=t$ for $i\in B_t$ and, given whatever $\sum _{i\in B_t}{w}_i$, the term $\sum _{i\in B_t}{w}_i^2$ is minimised at $w_i=w_t$ for $i\in B_t$.

As the first validity condition for $\hat{Y}$, suppose there exists a vector ${\text{β}}_{K\times 1}$, such that (18) $\sum _{i\in U_t}{ϵ}_i/N_t\to 0$(18) for each t-value, as $N\to \mathrm{\infty }$, where ${ϵ}_i=y_i-t_i^{\mathrm{\top }}\text{β}$, and $N_t$ is the population size of $U_t=\{i;t_i=t,i\in U\}$. The condition (18(18) $\sum _{i\in U_t}{ϵ}_i/N_t\to 0$(18) ) is analogous to the SP${}_x$ assumption (6(6) $E(y_i|x_i,i\in U)=\mu \left(x_i\right)=x_i^{\mathrm{\top }}\text{β},$(6) ), where the covariate $x_i$ is replaced by $t_i$ here. Moreover, it relaxes the model (6(6) $E(y_i|x_i,i\in U)=\mu \left(x_i\right)=x_i^{\mathrm{\top }}\text{β},$(6) ) of the conditional mean, allowing for potential heterogeneous mean similar to (15(15) $\mu \left(x\right)=\sum _{i\in U_x}{\mu }_i/N_x$(15) ). Now that $\sum _{i\in B}{w}_i{t}_i=T$, we have $\begin{array}{rl}\hat{Y}-Y& =\sum _{i\in B}{w}_i(t_i^{\mathrm{\top }}\text{β}+{ϵ}_i)-\sum _{i\in U}{t}_i^{\mathrm{\top }}(\text{β}+{ϵ}_i)\\ & =\sum _{i\in B}{w}_i{ϵ}_i-\sum _{i\in U}{ϵ}_i.\end{array}$ Given (18(18) $\sum _{i\in U_t}{ϵ}_i/N_t\to 0$(18) ), $\sum _{i\in U}{ϵ}_i/N\to 0$ as $N\to \mathrm{\infty }$. Moreover, we have $\begin{array}{rl}& \frac{1}{N}\sum _{i\in B}{w}_i{ϵ}_i=\sum _t\frac{w_t}{N}\sum _{i\in U_t}{\delta }_i{ϵ}_i\\ & =\sum _t{w}_t\frac{N_t}{N}\left({\mathrm{C}\mathrm{o}\mathrm{v}}_{N_t}({\delta }_i,{ϵ}_i)+\left(\frac{1}{N_t}\sum _{i\in U_t}{\delta }_i\right)\right.\\ & \left.\left(\frac{1}{N_t}\sum _{i\in U_t}{ϵ}_i\right)\right)\to 0\end{array}$ as $N\to \mathrm{\infty }$, given (19) $\begin{array}{rl}& {\mathrm{C}\mathrm{o}\mathrm{v}}_{N_t}({\delta }_i,{ϵ}_i)\to 0\\ & E_{N_t}\left({\delta }_i\right)=\sum _{i\in U_t}{\delta }_i/N_t\to p_t>0\end{array}$(19) for any given t, which is an adaption of the NPA non-informativeness assumption (17(17) $\begin{array}{rl}\underset{N\to \mathrm{\infty }}{lim}{\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)=0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ \underset{N\to \mathrm{\infty }}{lim}{E}_N\left({\delta }_i\right)=p>0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}$(17) ) to the present setting. It follows that the two assumptions (18(18) $\sum _{i\in U_t}{ϵ}_i/N_t\to 0$(18) ) and (19(19) $\begin{array}{rl}& {\mathrm{C}\mathrm{o}\mathrm{v}}_{N_t}({\delta }_i,{ϵ}_i)\to 0\\ & E_{N_t}\left({\delta }_i\right)=\sum _{i\in U_t}{\delta }_i/N_t\to p_t>0\end{array}$(19) ) are the key validity conditions for the calibration estimator to be consistent for Y.

For variance estimation, suppose again independent Bernoulli distribution of ${\delta }_i$ with probability $p_i$, where $0\le p_i\le 1$. An approximate variance estimator for the calibration estimator $\hat{Y}$ can then be given as $\hat{V}\left(\hat{Y}\right)=\sum _t({\hat{p}}_t^{-1}-1){\hat{p}}_t^{-1}\sum _{i\in B_t}(y_i-t_i^{\mathrm{\top }}\widehat{\text{β}}{)}^2,$ where ${\hat{p}}_t=n_{tB}/N_t$, and $\widehat{\text{β}}=\left(\sum _{i\in B}{w}_i{t}_i{t}_i^{\mathrm{\top }}{)}^{-1}\right(\sum _{i\in B}{w}_i{t}_i{y}_i)$.

3.4. Validation of non-informative B-sample selection

Of the validity conditions discussed above, the critical assumption is non-informative B-sample selection, which can be stated in various forms. For instance, given the non-informativeness assumption (17(17) $\begin{array}{rl}\underset{N\to \mathrm{\infty }}{lim}{\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)=0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ \underset{N\to \mathrm{\infty }}{lim}{E}_N\left({\delta }_i\right)=p>0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}$(17) ), an additional assumption like (18(18) $\sum _{i\in U_t}{ϵ}_i/N_t\to 0$(18) ) can in principle to validated empirically. However, the non-informativeness condition may not hold exactly, and it is generally impossible to verify only based on the data used for the estimation. Below we discuss the issue in more detail.

Consider first the propensity model $p_i=p(x_i;\eta )$ under the QR approach. Suppose known $x_U$ to avoid additional complications otherwise, the census score equation is $\sum _x\frac{\mathrm{\partial }p(x;\eta )}{\mathrm{\partial }\eta }\left[\frac{n_{xB}}{p(x;\eta )}-\frac{N_x-n_{xB}}{1-p(x;\eta )}\right]=0,$ which is always satisfied by $p(x;\hat{\eta })=n_{xB}/N_x$, i.e., the saturated model. Insofar as a non-saturated model of $p(x_i;\eta )$ does not fit perfectly to the data, one can always attribute its cause to the non-saturated functional form of $p(x_i;\eta )$, instead of rejecting the assumption that the set of covariates $x_i$ ‘fully govern the sampling mechanism’. In this sense, the validity of the latter assumption cannot be refuted empirically.

Next, for the SP approach, where both ${\delta }_i$ and $y_i$ are treated as random, assume the B-sample inclusion probability $p_i$ depends on $x_i$, where $x_i$ is known for $i\in U$ to avoid extra complications. For goodness-of-fit checks, let $z_i$ be a known covariate, which is distinct from $x_i$. We have $\begin{array}{rl}E\left(z_B\right)& =\sum _{i\in U}{p}_i{z}_i=\sum _{x}p(x;\eta )\sum _{i\in U_x}{z}_i\\ & =\sum _{x}p(x;\eta )N_z{\overline{Z}}_x,\\ Z& =E\left(\sum _{i\in U}{\delta }_i{z}_i/p_i\right)=E\left[\sum _x{n}_{xB}{\overline{z}}_{xB}/p(x;\eta )\right],\end{array}$ where ${\overline{Z}}_x=\sum _{i\in U_x}{z}_i/N_x$ and ${\overline{z}}_{xB}=\sum _{i\in B_x}{z}_i/n_{xB}$. The two observed checks are $\begin{array}{rl}& \left\{\begin{array}{l}{z}_B\equiv \sum _x{n}_{xB}{\overline{z}}_{xB}=\sum _x{\hat{p}}_x{N}_x{\overline{Z}}_x\\ Z=\sum _x{n}_{xB}{\overline{z}}_{xB}/{\hat{p}}_x\end{array}\right.\\ & \stackrel{\mathrm{i}\mathrm{f}{z}_i\equiv 1}{\Rightarrow }\left\{\begin{array}{l}\sum _{i\in U}{\hat{p}}_i=n_B\\ \sum _{i\in B}1/{\hat{p}}_i=N.\end{array}\right.\end{array}$ Setting ${\hat{p}}_x=n_{xB}/N_x$, which fits the assumption $p_i=p(x_i;\lambda )$, both the two checks are satisfied given ${\overline{Z}}_x={\overline{z}}_{xB}$, i.e., the B-sample expansion estimate of $Z_x$ is perfect for all x. This would suggest that the NPA assumption (17(17) $\begin{array}{rl}\underset{N\to \mathrm{\infty }}{lim}{\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)=0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ \underset{N\to \mathrm{\infty }}{lim}{E}_N\left({\delta }_i\right)=p>0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}$(17) ) holds for $z_i$ given $x_i$ and may be considered to support the plausibility of the NPA assumption (17(17) $\begin{array}{rl}\underset{N\to \mathrm{\infty }}{lim}{\mathrm{C}\mathrm{o}\mathrm{v}}_N({\delta }_i,y_i)=0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ \underset{N\to \mathrm{\infty }}{lim}{E}_N\left({\delta }_i\right)=p>0,& \mathrm{i}.\mathrm{e}.,\mathrm{n}\mathrm{o}\mathrm{n}\text{-}\mathrm{n}\mathrm{e}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{B}\text{-}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}$(17) ) for $y_i$ given $x_i$, provided $z_i$ is known to be correlated with $y_i$, but not otherwise. However, in situations where such a covariate $z_i$ is available, it seems natural that it should be used in the estimation of Y to start with. The two checks amount then to the case of $z_i\equiv 1$ and are satisfied trivially by setting ${\hat{p}}_x=n_{xB}/N_x$. Thus one is faced with a dilemma, where building the best model for estimation would at the same time reduces the ability to verify it.

4. Using additional probability sample of outcomes

So far we have only considered the situations, where the outcome values of interest are only observed in the non-probability sample B. Obviously, the situation changes completely, given in addition a probability sample of outcomes. Below we discuss shortly two different approaches to inference in the absence of any relevant covariates. The ideas remain the same in situations with additional covariates.

The first approach aims at consistent estimation combing the two samples, as for example, discussed in Tam and Kim (2018a, 2018b), where the probability sample is taken from the whole population and overlaps with the B-sample. These authors also discussed additional issues such as measurement errors or nonresponse. Here we discuss the situation where the probability sample is taken from the B-sample complement population. Given the non-probability sample observations $y_B$, one may treat $(B,y_B)$ as fixed, and select a second supplementary sample from the rest of the population, denoted by $S\subset U\setminus B$. Given the S-sample observations of the outcome, denoted by $y_S$, it is straightforward to obtain a test for $H_0:\overline{Y}={\overline{y}}_B$ vs. $H_1:\overline{Y}\ne {\overline{y}}_B$, given as $D=({\overline{y}}_B-\widehat{{\overline{Y}}_B^c}{)}^2/\hat{V}(\widehat{{\overline{Y}}_B^c})\sim {\chi }_1^2,$ where $\widehat{{\overline{Y}}_B^c}$ is an S-sample estimator of the population mean outside of the B-sample, i.e., ${\overline{Y}}_B^c=\sum _{i\in U\setminus B}{y}_i/(N-n_B)$ and $\hat{V}\left(\widehat{{\overline{Y}}_B^c}\right)$ is the associated variance estimator. If $H_0$ is not rejected, then there is the possibility of using ${\overline{y}}_B$ as an estimate on its own, without regular concurrent surveys in future. This would achieve the greatest cost savings. To this end, one may consider S as a particular form of audit sampling, whose aim is to validate the big-data estimate ${\overline{y}}_B$ and to provide a meaningful accuracy measure that can accommodate its potential bias. Zhang (2019) develops an approach to audit sampling inference for big data statistics.

Let $W_B=n_B/N$. A consistent estimator of $\overline{Y}$ using both samples is given by ${\hat{\overline{Y}}}_S=W_B{\overline{y}}_B+(1-W_B){\overline{y}}_w\mathrm{a}\mathrm{n}\mathrm{d}{\overline{y}}_w=\frac{\sum _{i\in S}{y}_i/{\pi }_i}{\sum _{i\in S}1/{\pi }_i},$ where ${\pi }_i$ is the S-sample inclusion probability, and the validity of ${\hat{\overline{Y}}}_S$ now derives from probability sampling of S, regardless of how the B-sample is generated. The relative efficiency (RE) against the setting without the B-sample can be given by $\text{RE}=[(1-W_B{)}^{2}V({\overline{y}}_w\left)]/V\right(\widehat{{\overline{Y}}^{\prime }}),$ where $\widehat{{\overline{Y}}^{\prime }}$ is a hypothetical probability sample from the whole population U, which has the same sample size and the same sampling design as S. One may refer to this as the split-population approach to inference, which is an age-old idea for combining survey sampling with administrative data. The efficiency gain would be substantial provided the B-sample is large. In fact, the larger the B-sample, the greater is the efficiency gain.

Under the second approach to inference, consider a composite estimator given by ${\hat{\overline{Y}}}_C=\gamma {\overline{y}}_B+(1-\gamma ){\overline{y}}_w,$ where γ is the composition weight, for $W_B\le \gamma \le 1$. Notice that when $\gamma =W_B$, the composite estimator is just the split-population estimator ${\hat{\overline{Y}}}_S$ above, which is consistent for $\overline{Y}$. As γ increases from $W_B$ towards one, one risks introducing greater bias, insofar as ${\overline{y}}_B\ne \overline{Y}$. However, the composite estimator may yield a smaller mean squared error (MSE) of estimation, provided this is desirable. One is then essentially trading the increasing bias $(\gamma -W_B)({\overline{y}}_B-{\overline{Y}}_B^c)$ against the decreasing stand error $(1-\gamma )\text{SE}\left({\overline{y}}_w\right)$, as γ increases. The composite estimator that achieves the minimum MSE is given by $\gamma =\frac{V\left({\overline{y}}_w\right)+W_B({\overline{y}}_B-{\overline{Y}}_B^c{)}^2}{V\left({\overline{y}}_w\right)+({\overline{y}}_B-{\overline{Y}}_B^c{)}^2}.$ Estimating ${\overline{Y}}_B^c$ by ${\overline{y}}_w$ in application, one can use $\hat{\gamma }=min(W_B+(1-W_B\left)\hat{V}\right({\overline{y}}_w\left)/\right({\overline{y}}_B-{\overline{y}}_w{)}^2,1).$ The validity of the composite approach derives from probability sampling of S, regardless of how the B-sample is generated. Again, the bigger the B-sample, the better it is.

5. Summary

All the estimators from non-probability sample observations reviewed in Section 2 are model based, whether the modelling is carried out under the SP or QR approach. Two features regarding the model covariate $x_i$, for $i\in U$, are worth recapitulating:

• compared to the situation with known $x_U$, making use of an additional probability sample $x_S$ entails a loss of efficiency, as can be expected;

• the availability of an additional probability sample without the outcome variable is not a principal advantage, since it does not simplify the validity conditions compared to the situation where $x_U$ is known, but it does resolve the practical difficulty when $x_U$ is unavailable yet some functions of $x_U$ are needed for descriptive inference.

The situation changes completely, given in addition a probability sample of outcomes. The probability sample then enables valid descriptive inference in combination with the non-probability probability sample. Depending on the circumstances, the probability sample can either be selected from the whole population, or just the rest population outside the non-probability sample.

Finally, in situations where the non-probability sample is large, the cost savings will be the greatest if it can replace regular survey sampling altogether. Use of an additional probability audit sample is needed to validate the non-probability sample estimate, in spite of possible failure of its underlying model assumptions, and to provide a meaningful accuracy measure that can accommodate its potential bias.

Disclosure statement

No potential conflict of interest was reported by the author.

ORCID

Li-Chun Zhang  http://orcid.org/0000-0002-3944-9484

Additional information

Notes on contributors

Li-Chun Zhang

Li-Chun Zhang is Professor of Social Statistics at University of Southampton, Senior Researcher at Statistics Norway, and Professor of Official Statistics at University of Oslo. His research interests include finite population sampling design and coordination, graph sampling, sample survey estimation, non-response, measurement errors, small area estimation, index number calculations, editing and imputation, register-based statistics, analysis of integrated data, statistical matching, record linkage, population size estimation.