Design-Based Predictive Inference

Abstract

Design-based inference from probability samples is valid by construction for target parameters that are descriptive summaries of finite populations. We develop a novel approach of design-based predictive inference for finite populations, where the individual-level predictor is learned from a probability sample using any models or algorithms for incorporating the relevant auxiliary information, and the uncertainty of estimation is evaluated with respect to the known probability design while the outcome and auxiliary values for modeling are treated as constants. Unlike the existing theory of design-based model-assisted estimation for finite populations, design-based predictive inference is as well suited for individual-level prediction in addition to producing population-level estimates.

1. Introduction

Throughout the twentieth century, design-based inference from finite-population probability sampling has been established as the standard approach to official statistics; see Hansen (1987), Smith (1994), Kalton (2002), Rao (2005, 2011), Beaumont and Haziza (2022) for reviews and appraisals. In this context the target parameters for estimation are descriptive, observable summaries of a given finite population, such as the population total, mean or quantiles of some specific values associated with the given population units, and the inference is characterized as descriptive or predictive (Smith, 1983; Geisser, 1993), in contrast to analytic inference of theoretical, unobservable targets such as the life expectancy (of a hypothetical cohort of individuals) or a parametric model that can be used to understand the given population.

By design-based inference from probability samples, the uncertainty of estimation is evaluated with respect to hypothetically repeated sampling from the same finite population, while all the other values involved are treated as constants associated with the given population. Design-based inference is valid by construction because it is based on the known sampling design, “whatever the unknown properties of the population” (Neyman 1934). In contrast, by model-based inference, the uncertainty of estimation is evaluated with respect to an assumed statistical model of the observations, while the available sample is typically treated as fixed; see for example, Valliant et al. (2000). Although models are necessary for analytic targets or if the available observations are not obtained by probability sampling, model-based inference may be invalid to the extent the assumed model is misspecified in respects that matter to the task at hand.

Design-based inference can be made more efficient by using auxiliary information in addition to the sampling design. This can largely recover the ‘loss of efficiency’ compared to model-based inference that uses the same auxiliary information optimally under the assumed model. For instance, calibration estimation (Deville and Särndal 1992) is a general approach that makes adjustments to the design weights with respect to the known auxiliary population totals. Or, empirical likelihood methods can yield confidence intervals that have better properties than normal approximation based on the central limit theorem (Hartley and Rao, 1968; Rao and Wu, 2010; Berger and De La Riva Torres, 2016). More relevant to our development is the model-assisted approach where an assisting model is explicitly formulated but inference remains design-based, whether or not the adopted estimator has optimal properties with respect to the assisting model. One can use linear models (Särndal et al. 1992), generalized linear models (Wu and Sitter 2001), or many other models under a unified “construction recipe” as reviewed by Breidt and Opsomer (2017).

To justify any model-assisted estimator that is not design-unbiased, it is common to seek a proof that it can be design consistent asymptotically for a hypothetical sequence of populations of increasing sizes. As Smith (1994) points out, this “asymptotic notion of consistency” is not immediately applicable to the given population as “a real entity”. In contrast, for a given population and sampling method, if $t(1),...,t(k)$ are unbiased estimators of the population totals $T(1),...,T(k)$ , then $g(t(1),...,t(k))$ is called “consistent” for $g(T(1),...,T(k))$ by Fisher (1956), in that replacing $t(j)$ by $T(j)$ gives the true target population parameter. Similarly, an interval estimator of “a collective character … of a population” is called “consistent” by Neyman (1934), if it achieves the designated level of coverage given the finite population and sampling method.

We emphasize that asymptotic consistency of a point estimator or an interval estimator is unnecessary, if the estimator is Neyman-Fisher consistent in the sense Neyman and Fisher have used the term “consistent” for finite-population inference. This will be our perspective to design-based predictive inference in this paper, which may also be called fully design-based inference in contrast to asymptotically design-consistent inference (that has been traditionally more common for model-assisted finite population estimation).

Now, we notice that design-unbiased ratio or linear regression estimators for population totals have been proposed by Hartley and Ross (1954) and Mickey (1959), which are finite-population Neyman-Fisher consistent. More recently, Sanguiao-Sande and Zhang (2021) developed a design-unbiased approach, called subsampling Rao-Blackwellisation (SRB), which allows for any assisting Machine Learning (ML) models or algorithms that have become increasingly common. The SRB approach combines three classic ideas in statistical inference, (i) model-assisted estimation for survey sampling, (ii) cross-validation for error estimation by ML methods, and (iii) the Rao-Blackwell Theorem (Rao, 1945; Blackwell, 1947) for efficiency improvement.

In this paper, we extend the SRB approach to a larger class of population estimators, which are commonly referred to as the prediction estimators, as well as the associated individual-level predictors for the out-of-sample units. Notice that, traditionally, due to the lack of a design-based prediction theory, individual outcomes must be treated as random variables for model-based prediction and the term predictor is common in this context. Although from a design-based inference perspective the term prediction estimator would seem more appropriate also at the unit level, we shall keep the term predictor at the individual level for convenience and familiarity reasons. Notice also that the two terms “unit” and “individual” are used interchangeably in this paper, such as in “statistical unit” or “individual prediction.”

It may be helpful to make some remarks immediately regarding the nature of our inference approach and its advantages compared to the more familiar model-assisted inference.

First, we consider predictive inference by definition, where the sample-based prediction estimator (using any given ML model or algorithm) aims at some out-of-sample quantity that varies with the sample, the property of which is evaluated only with respect to repeated sampling from the given population. This design-based predictive inference outlook differs from model-assisted estimation that is aimed at fixed population parameters (such as totals or means).

Next, we develop Neyman-Fisher consistent uncertainty estimators, which accommodate any given assisting model (or algorithm) and apply generally to sampling from finite populations. In contrast, asymptotically design-consistent inference may not hold in a given setting of finite population sampling, but still requires tailored asymptotic arguments to be developed for different nonparametric assisting models, such as random forest or support vector machine, which will remain a challenge as new models or algorithms emerge.

Finally, while our approach is model-assisted in the sense that the sampling design remains the inference basis despite the model introduced, it provides as well a design-theoretical basis for individual level prediction (estimation). This is another important difference to standard model-assisted estimation that is only applicable to population parameter estimation.

1.1. Prediction Estimator

Denote by $U=\{1,....,N\}$ a given finite population that is of size $N$ . Let $y_U=\{y_i:i\in U\}$ be the associated values of interest. Denote by $x_U=\{x_i:i\in U\}$ the collection of feature vectors, where $x_i$ is the vector associated with each unit $i\in U$ . Given any sample of units from $U$ , denoted by $s\subset U$ , let $\mu (x,s)$ be a predictor for any out-of-sample unit, say, $j\in R=U\backslash s$ , whose feature vector takes value $x$ , that is, $x_j=x$ . Notice that any function of $\{y_i:i\in U\backslash s\}$ is a random quantity that varies with the sample $s$ , just like $\mu (x,s)$ itself as the notation emphasizes. When the target parameter is the population total $Y=\sum _{i\in U}{y}_i$ , the prediction estimator of $Y$ is given as

$$\hat{Y}=\sum _{i\in s}{y}_i+\sum _{j\in U\backslash s}\mu (x_j,s)$$

(1)

We shall consider $\mu (x,s)$ as the associated individual-level predictor of $y$ for any unit with given features $x$ , where $y$ is treated as a constant just like $x$ .

Notice that individual level features are required to compute the prediction estimator, Equation (1), except for example, in the case of using linear models. However, this is a requirement common to all individual-level prediction models, regardless the inference framework. There may be situations where such information is unavailable, which would limit one’s choice of models. But the ability to utilize individual level covariates is certainly not a limitation of our approach not least because, as we shall explain, design-based individual-level predictive inference can provide a valid theoretical basis for producing census-like statistical data, which fills a gap in the current literature.

Note that many design-based estimators in survey sampling can as well be given as prediction estimators, Equation (1). For example, let $x_i={\pi }_{i}N/n$ for any $i\in U$ , given ${\pi }_i=\mathrm{Pr}(i\in s)$ and $n=|s|$ , such that the Horovitz-Thompson (HT) estimator (Horvitz and Thompson 1952) can be given in the form of Equation (1) with

$${\mu }_{\mathrm{HT}}(x_j,s)=x_j{\beta }_s+\frac{1}{N-n}\sum _{i\in s}(x_i{\beta }_s-y_i)$$

where ${\beta }_s=n^{-1}\sum _{i\in s}{y}_i/x_i$ . Other examples include the generalized regression estimator, model-calibrated linear estimator, as well as the SRB estimator of Sanguiao-Sande and Zhang (2021).

However, we are interested in design-unbiased inference of the prediction estimator, Equation (1) generally, including when $\mu$ is given by an arbitrary ML method without regard to the sampling design and may not have a wieldy expression, such as a random forest trained on $\{(y_i,x_i):i\in s\}$ by ready-made software.

The theory of design-based predictive inference for population totals and individuals will be developed and illustrated in Sections 2 and 3, respectively, an illustrative application to Structural Business Survey is given in Section 4, and some final remarks are given in Section 5. However, before we get into the details of the development, let us first motivate below what design-based predictive inference can do for official statistics.

1.2. Introduction to Total Estimation

Design-based predictive inference is clearly relevant to the perennial design versus model controversy in survey sampling, as traditionally the two main strands of approaches to finite population estimation.

In the design-based approach, estimators depend on the sampling design through the sample inclusion probabilities or other known sampling probabilities. Although auxiliary information in addition to the sampling design can be incorporated by various techniques, the validity and the associated uncertainty of the resulting estimator are still based on the given sampling design. In contrast, the model-based prediction approach, frequentist or Bayesian, depends on an assumed working model, which typically ignores the sampling design and treats the available sample as fixed when it comes to the assessment of the associated uncertainty.

Meanwhile, it is possible to evaluate any design-based estimator, such as the HT estimator or a generalized regression estimator, with respect to an assumed model, in which case it is common to conclude that design-based estimation is inefficient or lacks desirable conditional properties (e.g., Valliant et al. 2000). Conversely, a prediction estimator derived optimally under a working model can be evaluated with respect to the sampling design, in which case the danger of model misspecification is frequently noted (e.g., Hansen et al. 1983).

Design-based predictive inference takes the last analysis further, whereby one explicitly evaluates the design-based bias and variance of any given model-based prediction estimator, Equation (1). One can then compare the given model-based estimator to any other estimator, whether the latter is dependent on the design or a working model, and choose according to their design-based properties regardless how they are constructed. The merit of such an approach rests now on the fact that design-based uncertainty assessment is valid, which does not require the model underlying any given estimator to be correct.

To illustrate, under the linear model $E_M(y_i|x_i)=x^{\mathrm{T}}\beta$ and $V_M(y_i|x_i)={\sigma }^2$ where $E_M$ and $V_M$ denote expectation and variance under the model, the best linear unbiased predictor (BLUP) of a population total $Y=\sum _{i\in U}{y}_i$ is

$$\hat{Y}=X^{\mathrm{T}}b\mathrm{and}b=(\sum _{i\in s}{x}_i{x}_i^{\mathrm{T}}{)}^{-1}\sum _{i\in s}{x}_i{y}_i$$

where $X=\sum _{i\in U}{x}_i$ . Given ${\pi }_i=\mathrm{Pr}(i\in s)$ , we have

$$E_p(\hat{Y})=X^{\mathrm{T}}{E}_p(b)=Y-\sum _{i\in U}\{y_i-x_i^{\mathrm{T}}{E}_p(b)\}$$

where $E_p(b)\approx (\sum _{i\in U}{\pi }_i{x}_i{x}_i^{\mathrm{T}}{)}^{-1}\sum _{i\in U}{\pi }_i{x}_i{y}_i$ with respect to sampling, and

$$V_p(\hat{Y})=X^{\mathrm{T}}{V}_p(b|s)X\approx X^{\mathrm{T}}{E}_p\{V_M(b|s)\}X$$

with respect to sampling, where

$$E_p\{V_M(b|s)\}=E_p\left\{{\sigma }^2{\left(\sum _{i\in s}{x}_i{x}_{i^{\top }}\right)}^{-1}\right\}\approx \left(\frac{1}{N}\sum _{i\in U}{\{y_i-x_i^{\mathrm{T}}{E}_p(b)\}}^2\right){\left(\sum _{i\in U}{\pi }_i{x}_i{x}_i^{\mathrm{T}}\right)}^{-1}.$$

Both $E_p(\hat{Y})$ and $V_p(\hat{Y})$ can then be estimated from $s$ and compared to, say, the approximate variance of a generalized regression estimator of $Y$ .

Thus, a chief advantage of adopting design-based predictive inference to population total estimation based on probability sampling is to circumvent the design versus model controversy, by providing a valid common ground for uncertainty assessment. A theory applicable to the class of prediction estimators, Equation (1), which will be developed in this paper, would allow one to use any assisting ML models or algorithms that can often make more efficient use of auxiliary information than the standard design-based calibration estimation or model-assisted estimation methods.

1.3. Introduction to Individual Estimation

Individual estimation requires the most extremely disaggregated results. It can be useful for constructing statistical registers or census-like statistical data as the basis for descriptive official statistics. However, there has never been a design-based theory for estimation at the individual level.

For instance, having taken a simple random sample of all but two units in a given population, the traditional design-based estimation theory would only allow one to make inference about the total (or mean) of the two out-of-sample units, but not each on its own, no matter how large the sample is or how much auxiliary information one has in addition. This is clearly unsatisfactory, which requires extension of the design-based inference theory.

To illustrate the conceptual issue at hand, suppose on observing $\{y_i:i\in s\}$ in a subset $s$ of the population $U$ , one would like to predict the $y$ -value for each unit out of $s$ by $\mu (s)=\sum _{i\in s}{y}_i/n$ , where $n$ is the number of units in $s$ . How can one infer about the loss $D_s=\sum _{j\in U\backslash s}\{\mu (s)-y_j{\}}^2$ that is unobserved?

One possibility is to assume a model. For instance, under the model that $y_i$ is independent and identically distributed (IID) for any $i\in U$ , we have

$$E_M(D_s|s)=(N-n)(1+n^{-1}){\sigma }^2$$

with respect to the IID model conditional on the given subset $s$ , where ${\sigma }^2$ is the variance of $y_i$ under the model and $N$ is the number of units in $U$ .

However, we notice that a fundamentally different, design-based approach would in fact be possible if $s$ is selected from $U$ by a known sampling design, denoted by $s~p(s)$ , where $\sum _{s\in \Omega }p(s)=1$ and $\Omega$ contains all the possible samples from $U$ . For instance, suppose $s$ is selected from $U$ by simple random sampling without replacement (SRSWOR), where $p(s)=1/\left(\begin{array}{c}N\hfill \\ \multicolumn{1}{c}{n}\end{array}\right)$ , such that

$$E_p(D_s)=(N-n)(1+n^{-1})S_y^2$$

with respect to $p(s)$ , where $S_y^2=\sum _{i\in U}(y_i-\overline{Y}{)}^2/(N-1)$ and $\overline{Y}=\sum _{i\in U}{y}_i/N$ .

Since $s_y^2=\sum _{i\in s}\{y_i-\mu (s){\}}^2/(n-1)$ is both unbiased for ${\sigma }^2$ under the IID model and unbiased for $S_y^2$ under SRSWOR, numerically one would obtain the same estimate of the expected loss, even though they have completely different interpretations. While an assumed model would be necessary for evaluating $E_M(D_s|s)$ if the selection mechanism of $s$ is unknown, it could be invalid if the observed data distribution actually differs to that of the unobserved ones. While the design-based loss $E_p(D_s)$ here requires one to plan and implement the SRSWOR design, it is necessarily valid because $p(s)$ is known.

We shall develop a general design-based theory for the out-of-sample loss, provided the observations are obtained by probability sampling. This would yield valid inference of the associated risk with respect to the given sampling design, where all the outcomes $y_U$ and features $x_U$ are treated as constants.

2. Total Prediction Estimator

Consider the prediction estimator, Equation (1), where $\mu (x,s)$ can be obtained by any model or algorithm fitted to the full sample $s$ . Now, hypothetically speaking, it is clear that design-unbiased estimation of $\hat{Y}-Y=\sum _{j\in R}\mu (x_j,s)-y_j$ , or some function of it, would be possible given an additional probability sample $r$ selected from $U\backslash s$ , because one can then observe the error $e_j=\mu (x_j,s)-y_j$ for any $j\in r$ . In the absence of extra observations $\{y_j:j\in r\}$ , valid design-based inference requires creating observed errors within the sample $s$ .

Denote by $s_1\cup s_2=s$ and $s_1\cap s_2=\varnothing$ a training-test sample split, where $s_1$ is selected by a subsampling design, denoted by $q(s_1|s)$ . For example, $s_1$ of size $n_1$ can be randomly sampled from $s$ with or without replacement. Or, in $T$ -fold cross-validation, $s$ is first randomly partitioned into $T$ clusters and then each cluster is selected as $s_2$ one by one systematically, yielding $s_1=s\backslash s_2$ accordingly. Denote by $\mu (x,s_1)$ the predictor obtained from the subsample $s_1$ , in the same way as $\mu (x,s)$ from $s$ . Its error $\mu (x_j,s_1)-y_j$ can be observed for any $j\in s_2$ .

As in Sanguiao-Sande and Zhang (2021), we shall refer to the sampling design that yields $(s_1,s)$ as the $\mathrm{pq}$ -design, denoted by

$$f(s_1,s)=q(s_1|s)p(s)=f(s|s_1)f(s_1)$$

(2)

where the last product indicates that, conditional on the training set $s_1$ , one can view the test set $s_2$ as a probability sample from $U\backslash s_1$ , according to which $s$ can vary under the $\mathrm{pq}$ -design. In particular, for any $i\in U$ , let

$${\pi }_{2i}=\mathrm{Pr}(i\in s_2|s_1)=\sum _{s\in i,i\notin s_1}f(s|s_1)$$

(3)

be its conditional $s_2$ -inclusion probability given $s_1$ under the $\mathrm{pq}$ -design.

The theory of design-based predictive inference we develop below, including both Theorems 1 and 2, apply generally provided any well-defined $\mathrm{pq}$ -design, Equation (2). However, the typical $q$ -designs mentioned above may need to be modified in practice, in order to accommodate the additional complex sampling features that may exist. For instance, given a stratified $p$ -design of $s$ , it may be natural to subsample $s_1$ within each stratum as well. Or, given a multistage $p$ -design, the $q$ -design must involve subsampling of the selected primary sampling units (PSUs), instead of only subsampling elements within all the selected PSUs, such that conditional sampling by $f(s|s_1)$ covers the whole population.

Now, conditional on any given $s_1$ , the design-based bias and mean squared error (MSE) of the prediction estimator using $\mu (x,s_1)$ can be easily derived. The theory below explains how one can infer the design-based bias and MSE of the prediction estimator (1) that is trained using all the observations in $s$ , with respect to $s~p(s)$ , by appropriate averaging over the $q$ -design.

2.1. SRB Prediction Estimator

Given $s_1$ under the $\mathrm{pq}$ -design, the subsample-trained prediction estimator is

$${\hat{Y}}_1^{*}=\sum _{i\in s}{y}_i+\sum _{j\in U\backslash s}\mu (x_j,s_1)$$

whose total error for $Y$ is given by

$$B=\sum _{j\in U\backslash s}{e}_{1j}=B_1-B(s_2)$$

where $e_{1j}=\mu (x_j,s_1)-y_j$ for any $j\notin s_1$ , and

$$B_1=\sum _{j\in U\backslash s_1}{e}_{1j}\mathrm{and}B(s_2)=\sum _{j\in s_2}{e}_{1j}$$

Note that $B$ and $B(s_2)$ vary with $(s,s_2)$ conditional on $s_1$ while $B_1$ is fixed.

Applying Rao-Blackwellisation to ${\hat{Y}}_1^{*}$ yields a corresponding SRB prediction estimator in the form of Equation (1), which is given by

$${\hat{Y}}^{\mathrm{RB}}=\sum _{i\in s}{y}_i+\sum _{j\in U\backslash s}\overline{\mu }(x_j,s)$$

(4)

where

$$\overline{\mu }(x_j,s)=E_q\{\mu (x_j,s_1)|s\}$$

(5)

Notice that $\overline{\mu }(x,s)$ is a particular full-sample trained predictor, and its special notation $\overline{\mu }$ is introduced to distinguish it from any $\mu (x,s)$ that uses the same $\mu$ but is directly trained once on the full sample $s$ . Sanguiao-Sande and Zhang (2021) refer to the operation $E_q(·|s)$ as SRB, since it yields the conditional expectation over subsampling $s_1~q(s_1|s)$ , where the unordered $s$ is the minimal sufficient statistic with respect to the sampling distribution $p(s)$ .

It follows from the definition that the bias of ${\hat{Y}}^{\mathrm{RB}}$ is given by

$$E_p({\hat{Y}}^{\mathrm{RB}})-Y=E_{\mathrm{pq}}({\hat{Y}}_1^{*})-Y=E_{\mathrm{pq}}(B).$$

Given $s_1$ , a conditionally unbiased predictor of $B$ is given by

$$\hat{B}=\sum _{j\in s_2}({\pi }_{2j}^{-1}-1)e_{1j}$$

in the sense that, as $s$ varies conditional on $s_1$ , we have

$$E_s(\hat{B}|s_1)=B_1-E_s\{B(s_2)|s_1\}=E_s(B|s_1)$$

since $e_{1j}$ resulting from $\mu (x_j,s_1)$ trained on the subsample $s_1$ is fixed with respect to sampling of $s_2$ by $f(s|s_1)$ conditional on $s_1$ . Applying Rao-Blackwellisation to $\hat{B}$ yields then a more efficient estimator

$${\hat{B}}^{\mathrm{RB}}=E_q(\hat{B}|s)$$

(6)

as an unbiased estimator of the bias of ${\hat{Y}}^{\mathrm{RB}}$ with respect to $p(s)$ , since

$$E_p({\hat{B}}^{\mathrm{RB}})=E_p\{E_q(\hat{B}|s)\}=E_{s_1}\{E_s(\hat{B}|s_1)\}=E_{s_1}\{E_s(B|s_1)\}=E_{\mathrm{pq}}(B).$$

Notice that one needs at least $|s_2|=1$ to calculate $\hat{B}$ , in which case $\mu (x,s_1)$ is trained on $n-1$ units, that is, the so-called leave-one-out (LOO) predictor, whose difference to $\mu (x,s)$ is only due to one randomly selected unit.

Moreover, one can estimate unbiasedly the MSE of ${\hat{Y}}^{\mathrm{RB}}$ by the following result, the proof of which is given in Appendix.

• Theorem 1. For any given $\mu (·)$ , an unbiased estimator of the MSE of the SRB prediction estimator (4), over $s~p(s)$ , is given by

$${\mathrm{mse}}^{\mathrm{RB}}=E_q\{{\hat{B}}^2-{\hat{V}}_s(\hat{B}|s_1)+{\hat{V}}_s\{B(s_2)|s_1\}|s\}-V_q({\hat{Y}}_1^{*}|s)$$

(7)

where $\hat{B}=\sum _{j\in s_2}({\pi }_{2j}^{-1}-1)\{\mu (x_j,s_1)-y_j\}$ , and ${\hat{V}}_s(\hat{B}|s_1)$ is unbiased for

$$V_s(\hat{B}|s_1)=\sum _{i\notin s_1}\sum _{j\notin s_1}({\pi }_{2\mathrm{ij}}-{\pi }_{2i}{\pi }_{2j})(\frac{1}{{\pi }_{2i}}-1)(\frac{1}{{\pi }_{2j}}-1)e_{1i}{e}_{1j}$$

where ${\pi }_{2\mathrm{ij}}=\mathrm{Pr}(i,j\in s_2|s_1)$ , and ${\hat{V}}_s\{B(s_2)|s_1\}$ is unbiased for

$$V_s\{B(s_2)|s_1\}=\sum _{i\notin s_1}\sum _{j\notin s_1}({\pi }_{2\mathrm{ij}}-{\pi }_{2i}{\pi }_{2j})e_{1i}{e}_{1j}.$$

Notice that one needs at least $|s_2|=2$ to calculate the variance estimators in Equation (7), in which case $\mu (x,s_1)$ is trained on $n-2$ units, that is, the leave-two-out (LTO) predictor, which differs to $\mu (x,s)$ only due to two randomly selected units.

2.2. Discussion

In a recent discussion of cross-validation for prediction error estimation under the IID model of $(y_i,x_i)$ , Bates et al. (2024, Theorem 3) have shown that, given a sample $s$ of size $n$ , it is possible to estimate unbiasedly the MSE of $(K-1)$ -fold cross-validation on a reduced sample of size $n(K-1)/K$ , but not for the $K$ -fold cross-validation on the sample of size $n$ which one actually does. Unbiased MSE estimation for the latter is generally difficult if $\mu (x,s)$ does not have a wieldy expression, because by definition one cannot observe the error of $\mu (x,s)$ without extra out-of-sample observations.

We have obtained design-unbiased MSE estimator, Equation (7), for the SRB prediction estimator ${\hat{Y}}^{\mathrm{RB}}$ that uses $\overline{\mu }(x,s)$ , but not for $\hat{Y}$ using the more familiar $\mu (x,s)$ that is trained once on the full sample. Design-unbiased estimation of $\mathrm{MSE}(\hat{Y})$ faces the same general difficulty as model-based MSE estimaton. Indeed, similarly to reduced-sample cross-validation, we could easily estimate unbiasedly the design-based MSE of the hypothetical prediction estimator

$${\hat{Y}}_1=\sum _{i\in s_1}{y}_i+\sum _{j\in U\backslash s_1}\mu (x_j,s_1)$$

as if $s_1$ was the sample (instead of the actual $s$ ) such that $\mu (x,s_1)$ was the full-sample once-trained predictor. Analogously to Theorem 1, we would have $E_q\{{\hat{B}}_1|s\}$ as an unbiased estimator of the design-based bias of ${\hat{Y}}_1$ , where

$${\hat{B}}_1=\sum _{i\in s_2}{\pi }_{2i}^{-1}\{\mu (x_i,s_1)-y_i\}$$

and $E_q\{{\hat{B}}_1^2-{\hat{V}}_s({\hat{B}}_1)|s\}$ unbiasedly for $\mathrm{MSE}({\hat{Y}}_1)$ , where ${\hat{V}}_s({\hat{B}}_1)$ is unbiased for

$$V_s({\hat{B}}_1)=\sum _{i\notin s_1}\sum _{j\notin s_1}(\frac{{\pi }_{2\mathrm{ij}}}{{\pi }_{2i}{\pi }_{2j}}-1)\{\mu (x_i,s_1)-y_i\}\{\mu (x_j,s_1)-y_j\}$$

Next, as will be illustrated later, we note that $\mathrm{MSE}({\hat{Y}}^{\mathrm{RB}})$ using $\overline{\mu }(x,s)$ can provide a good approximation to $\mathrm{MSE}(\hat{Y})$ using $\mu (x,s)$ , where ${\hat{Y}}^{\mathrm{RB}}$ and $\hat{Y}$ should be close to each other now that both of them are trained on the full sample, as long as $\mu$ is “stable” in a suitable sense. Specifically, we have

$$\hat{Y}-Y\equiv E_q(\hat{Y}-Y|s)=E_q\{\sum _{j\notin s}\mu (x_j,s)-\mu (x_j,s_1)|s\}+({\hat{Y}}^{\mathrm{RB}}-Y)$$

where the term $E_q\{·\}$ is of a lower order to ${\hat{Y}}^{\mathrm{RB}}-Y$ , if $\mu$ is $n_2$ -times $q$ -stable. Sanguiao-Sande and Zhang (2021) define the LOO-predictor $\mu$ to be $q$ -stable, if $\mu (x,s)-\mu (x,s_1){\to }^{P}0$ asymptotically, as $n,N\to \infty$ , while $n_2=1$ is held fixed. Here, $\mu$ is said to be $n_2$ -times $q$ -stable if the same holds given any $n_2\ge 1$ .

Finally, notice that, apart from familiarity and custom, there is no reason why one cannot adopt the SRB prediction estimator, Equation (4), for which MSE estimation is unbiased, instead of using the once-trained predictor $\mu (x,s)$ .

2.3. Notes on Implementation

2.3.1. Sampling Probability

By definition, Equation (3), ${\pi }_{2i}$ is the conditional $s_2$ -inclusion probability given $s_1$ , the calculation of which generally requires $f(s|s_1)$ that is derived from $q(s_1|s)p(s)$ . However, $p(s)$ is unknown for many unequal-probability sampling methods in practice, such as the cube method (Deville and Tillé 2004), although the inclusion probability ${\pi }_i=\mathrm{Pr}(i\in s)$ is always known.

One can use instead another sampling probability of the $\mathrm{pq}$ -design. For any $i\in U$ , let its conditional test-set inclusion probability given $i\notin s_1$ be

$${\varphi }_{2i}=\mathrm{Pr}(i\in s_2|i\notin s_1)=\frac{\mathrm{Pr}(i\in s_2,i\notin s_1)}{\mathrm{Pr}(i\notin s_1)}=\frac{{\pi }_i\{1-\mathrm{Pr}(i\in s_1|i\in s)\}}{1-{\pi }_i\mathrm{Pr}(i\in s_1|i\in s)}.$$

(8)

Given ${\pi }_i$ , the probability ${\varphi }_{2i}$ can be calculated as long as $\mathrm{Pr}(i\in s_1|i\in s)$ does not depend on $i$ under the subsampling design and can be specified regardless the realized $s$ , such as SRS of $s_1$ from $s$ with or without replacement, or $T$ -fold cross-validation. Since

$${\varphi }_{2i}=\frac{\sum _{s_1:i\notin s_1}\sum _{s:i\in s\cap i\notin s_1}f(s|s_1)f(s_1)}{\sum _{s_1:i\notin s_1}f(s_1)}=\frac{\sum _{s_1:i\notin s_1}{\pi }_{2i}f(s_1)}{\sum _{s_1:i\notin s_1}f(s_1)}=E_{s_1}\{{\pi }_{2i}|i\notin s_1\},$$

it is the conditional expectation of non-zero ${\pi }_{2i}$ , where ${\pi }_{2i}=0$ iff $i\in s_1$ .

To illustrate, take the special case of SRSWOR of $s$ from $U$ and SRSWOR of $s_1$ from $s$ , with sample sizes $n=|s|$ , $n_1=|s_1|$ , and $n_2=|s_2|=n-n_1$ . For any given $s_1$ and $i\notin s_1$ , we have exactly

$${\pi }_{2i}=\frac{\mathrm{Pr}(i\in s_2)f(s_1|i\in s_2)}{f(s_1)}=\frac{\frac{n_2}{N}{\left(\begin{array}{c}N-1\hfill \\ \multicolumn{1}{c}{n_1}\end{array}\right)}^{-1}}{{\left(\begin{array}{c}N\hfill \\ \multicolumn{1}{c}{n_1}\end{array}\right)}^{-1}}=\frac{n_2}{N-n_1}={\varphi }_{2i}.$$

Similarly, instead of ${\pi }_{2\mathrm{ij}}$ in Equation (7), one can use another conditional joint $s_2$ -inclusion probability

$${\varphi }_{2\mathrm{ij}}=\mathrm{Pr}(i,j\in s_2|i,j\notin s_1)=\frac{\mathrm{Pr}(i,j\in s_2|i,j\in s){\pi }_{\mathrm{ij}}}{\mathrm{Pr}(i,j\notin s_1)}$$

where $\mathrm{Pr}(i,j\in s_2|i,j\in s)$ is known for random subsampling of $s_1$ or $T$ -fold cross-validation, and $\mathrm{Pr}(i,j\notin s_1)=1-\mathrm{Pr}(i\in s_1)-\mathrm{Pr}(j\in s_1)+\mathrm{Pr}(i,j\in s_1)$ .

2.3.2. Monte Carlo SRB

By Equation (7), the exact ${\mathrm{mse}}^{\mathrm{RB}}$ for the LTO-SRB prediction estimator ${\hat{Y}}^{\mathrm{RB}}$ requires $\mathcal{C}(n,2)=\left(\begin{array}{c}n\hfill \\ \multicolumn{1}{c}{2}\end{array}\right)$ sample splits. If feasible, it would also provide a good estimator of $\mathrm{MSE}(\hat{Y})$ given any two-times $q$ -stable $\mu$ as explained earlier. The variance of ${\mathrm{mse}}^{\mathrm{RB}}$ should be comparable to that of the HT variance estimator; indeed, the former tends to be smaller than the latter if the model-assisted ${\hat{Y}}^{\mathrm{RB}}$ is more efficient than the HT-estimator ‘assisted’ only by $x_i={\pi }_{i}N/n$ (Section 1).

The exact SRB operation may be infeasible if $\mathcal{C}(n,2)$ is too large, in which case it needs to be replaced by Monte Carlo Rao-Blackwellisation (MC-RB) using $T$ sample splits, $T<<\mathcal{C}(n,2)$ . The corresponding MC-SRB prediction estimator, denoted by ${\tilde{Y}}^{\mathrm{RB}}$ , is given by replacing $\overline{\mu }(x,s)$ in Equation (4) by

$$\tilde{\mu }(x,s)=\frac{1}{T}\sum _{t=1}^T\mu (x,s_1^{(t)})$$

based on sample splits $(s_1^{(t)},s_2^{(t)})$ for $t=1,...,T$ . The MSE of ${\tilde{Y}}^{\mathrm{RB}}$ is then

$$\mathrm{MSE}({\tilde{Y}}^{\mathrm{RB}})=\mathrm{MSE}({\hat{Y}}^{\mathrm{RB}})+E_p\{V_q({\hat{Y}}_1^{*}|s)\}/T.$$

The unbiased exact-RB estimator of $\mathrm{MSE}({\tilde{Y}}^{\mathrm{RB}})$ follows from Equation (7), that is,

$${\mathrm{mse}}^{\mathrm{RB}}({\tilde{Y}}^{\mathrm{RB}})=E_q\{{\hat{B}}^2-{\hat{V}}_s(\hat{B}|s_1)+{\hat{V}}_s\{B(s_2)|s_1\}|s\}-\frac{T-1}{T}{V}_q({\hat{Y}}_1^{*}|s),$$

such that the unbiased MC-RB estimator of $\mathrm{MSE}({\tilde{Y}}^{\mathrm{RB}})$ is given by

$${\widetilde{\mathrm{mse}}}^{\mathrm{RB}}=\frac{1}{T}\sum _{t=1}^T({\hat{B}}^{2(t)}-{\hat{V}}_s(\hat{B}|s_1^{(t)})+{\hat{V}}_s\{B(s_2)|s_1^{(t)}\}-\{{\hat{Y}}_1^{*(t)}-{\tilde{Y}}^{\mathrm{RB}}{\}}^2),$$

(9)

where the index $t$ explicates the computation results for $t=1,...,T$ .

For the LTO ${\tilde{Y}}^{\mathrm{RB}}$ , the ratio $V_{\mathrm{pq}}({\widetilde{\mathrm{mse}}}^{\mathrm{RB}})/V_p({\mathrm{mse}}^{\mathrm{RB}})$ converges to 1 much slower than $\mathrm{MSE}({\tilde{Y}}^{\mathrm{RB}})/\mathrm{MSE}({\hat{Y}}^{\mathrm{RB}})$ , as $T$ increases. The inflation of $V_{\mathrm{pq}}({\widetilde{\mathrm{mse}}}^{\mathrm{RB}})$ is almost entirely due to approximating the $E_q\{·\}$ -term in Equation (7) by MC in Equation (9), wherein the terms are all estimated from $s_2^{(t)}$ of the size $n_2=2$ . While the MC variance of $V_{\mathrm{pq}}({\widetilde{\mathrm{mse}}}^{\mathrm{RB}})$ can be reduced by increasing either $T$ or $n_2$ , the reduction is considerably faster as $n_2$ increases. Note that the condition of $n_2$ -times $q$ -stability of $\mu$ is put under a greater pressure as $n_2$ increases, which might affect whether $\mathrm{MSE}({\tilde{Y}}^{\mathrm{RB}})$ would remain close to $\mathrm{MSE}(\hat{Y})$ , where $\hat{Y}$ is based on $\mu (x,s)$ .

Therefore, in practice, one can first explore the MC variance of ${\widetilde{\mathrm{mse}}}^{\mathrm{RB}}$ in relation to $n_2$ using a moderately large $T$ , in order to choose a value of $n_2$ that reduces the MC variance as much as possible without jeopardizing the $n_2$ -times stability condition, and then compute ${\widetilde{\mathrm{mse}}}^{\mathrm{RB}}$ with the chosen $n_2$ using $T$ that is as large as computationally practical.

2.4. Illustration

Let us illustrate with a simple example. Generate and fix a population of size $N=1000$ by $y_i={\beta }_1{x}_{1i}+{\beta }_2{x}_{2i}+{ϵ}_i$ with IID $x_{1i}~\mathrm{LogN}(1,1)$ , $x_{2i}~\mathrm{Poisson}(5)$ , and ${ϵ}_i~N(0,{\sigma }^2/4)$ , where ${\sigma }^2$ is the population variance of $x_{1i}$ . Let $s$ be given by SRSWOR with $n=100$ . Let SRSWOR be the subsampling $q$ -design of $s_1$ with size $n_1$ that is to be specified. Let the mis-specified predictor be $\mu (x,s)=x^{\mathrm{T}}\beta$ , with $x=(1,x_1)$ but omitting $x_2$ . The full-sample once-trained $\mu (x,s)$ and the SRB $\overline{\mu }(x,s)$ are approximately but not exactly equal to each other, where

$$\mu (x,s)=x^{\mathrm{T}}(\sum _{i\in s}{x}_i{x}_i^{\mathrm{T}}{)}^{-1}(\sum _{i\in s}{x}_i{y}_i)$$

$$\overline{\mu }(x,s)=E_q\{\mu (x,s_1)|s\}=x^{\mathrm{T}}{E}_q({(\sum _{i\in s}\mathbb{I}(i\in s_1)x_i{x}_i^{\mathrm{T}})}^{-1}(\sum _{i\in s}\mathbb{I}(i\in s_1)x_i{y}_i)|s)$$

Given a single realized sample $s$ , as in practice, Figure 1 illustrates the MC variance (left) and expectation (right) as $n_2$ increases from 2 toward $n$ under this setup, given $T={10}^3$ . The MC variance is seen to decrease sharply to a plateau as $n_2$ increases from 2, before it increases dramatically again as $n_2$ gets close to $n$ . Meanwhile, the MC expectation is quite stable, say, for $n_2\le 40$ . For example, setting $n_2$ to be 20 or 30 for the final computation of ${\widetilde{\mathrm{mse}}}^{\mathrm{RB}}$ using $T={10}^5$ , we would obtain $0.020$ or $0.014$ as the MC coefficient of variance (CV), which seem acceptable for practical purposes.

Figure 1.

MC variance (left) and expectation (right) of ${\widetilde{\mathrm{mse}}}^{\mathrm{RB}}$ given $n_2\ge 2$ .

Table 1 shows the results of simulating MSE estimation based on 250 independent samples, given $T={10}^3$ and $n_2=2,20,30$ . The MSE is simply the average squared error of either $\hat{Y}$ or ${\hat{Y}}^{\mathrm{RB}}$ over the 250 samples, and the relative efficiency (RE) is the ratio between either MSE and the variance of the HT-estimator. Notice that the three $\mathrm{MSE}(\hat{Y})$ here are all estimators of the same MSE, each using 250 independent samples, since $\hat{Y}$ depends only on $s$ .

Table 1.

MSE Estimation from 250 Samples, $T={10}^3$ , $\mu (x,s)$ for $\hat{Y}$ , and $\overline{\mu }(x,s)$ for ${\hat{Y}}^{\mathrm{RB}}$ (Training, Test) Set of Size $(n_1,n_2)$ , RE Against Variance of HT-Estimator.

$(n_1,n_2)$	$\mathrm{MSE}(\hat{Y})$	$\mathrm{RE}(\hat{Y})$	$\mathrm{MSE}({\hat{Y}}^{\mathrm{RB}})$	$\mathrm{RE}({\hat{Y}}^{\mathrm{RB}})$	$\mathrm{CV}({\widetilde{\mathrm{mse}}}^{\mathrm{RB}})$
(98, 2)	386,532.7	0.44	386,632.4	0.44	3.48
(80, 20)	363,613.9	0.41	363,441.5	0.41	0.31
(70, 30)	362,673.0	0.41	357,146.9	0.41	0.21

For $n_2$ up to 20 (or even 30), $\mathrm{MSE}({\hat{Y}}^{\mathrm{RB}})$ is practically equal to $\mathrm{MSE}(\hat{Y})$ . The CV of the MC-MSE estimator ${\widetilde{\mathrm{mse}}}^{\mathrm{RB}}$ is drastically reduced by setting $n_2$ to 20 or 30 instead of 2. In comparison, the CV of the exact-RB MSE estimator ${\mathrm{mse}}^{\mathrm{RB}}$ is 0.14 by simulation, whereas the CV of the HT variance estimator is 0.32. This confirms that setting $n_2$ to be 20 (or even 30) and using a larger but practical $T$ would work satisfactorily for MSE estimation in this setup.

In terms of the choice of estimator, we notice that the mis-specified predictor $\mu (x,s)$ yields a design-based MSE that is less than half of the variance of the HT-estimator, and the bias of $\hat{Y}$ or ${\hat{Y}}^{\mathrm{RB}}$ is a negligible part of the MSE here, the details of which are omitted to save space. The assessment is enabled by the design-based predictive inference theory developed above. Finally, as mentioned before, there is no reason why one cannot adopt ${\hat{Y}}^{\mathrm{RB}}$ using $\overline{\mu }(x,s)$ , for which MSE estimation is unbiased, instead of $\hat{Y}$ using $\mu (x,s)$ .

3. Individual Prediction Estimator

Consider the individual-level predictor $\mu (x,s)$ in the prediction estimator (1). Regardless how $\mu (x,s)$ is obtained from $\{(y_i,x_i):i\in s\}$ , using whichever model or algorithm, its total squared error (TSE) over $R=U\backslash s$ is given by

$$D(s;\mu )=\sum _{i\in R}{\{\mu (x_i,s)-y_i\}}^2.$$

For design-based individual-level predictive inference, we define the risk of $\mu$ to be the expectation of $D(s;\mu )$ over repeated sampling of $s~p(s)$ , while treating $y_U$ and $x_U$ as fixed, denoted by

$$\tau (\mu )=E_p\{D(s;\mu )\}.$$

(10)

We stress that only $s$ is random in Equation (10), that is, it is a design-based measure, regardless the model or algorithm by which $\mu (x,s)$ is constructed.

3.1. Risk of SRB Predictor

Under the same $\mathrm{pq}$ -design of $(s_1,s)$ as in Section 2, the error of the subsample-trained predictor $\mu (x,s_1)$ for any $i\notin s_1$ is given by

$$e_i(\mu ,s_1)=\mu (x_i,s_1)-y_i$$

which can be observed for any unit in the test set $s_2=s\backslash s_1$ . Let

$$D_R(s_1;\mu )=\sum _{i\in R}{e}_i(\mu ,s_1{)}^2$$

be the TSE of $\mu (x,s_1)$ over $R=U\backslash s$ . Let

$$A_2=\sum _{i\in s_2}{e}_i(\mu ,s_1{)}^2=\sum _{i\in U\backslash s_1}{e}_i(\mu ,s_1{)}^2-D_R(s_1;\mu ).$$

Given $s_1$ , both $A_2$ and $D_R(s_1;\mu )$ vary with $s_2=s\backslash s_1$ under the $\mathrm{pq}$ -design, but their sum $\sum _{i\in U\backslash s_1}{e}_i(\mu ,s_1{)}^2$ is fixed. The predictor

$${\hat{D}}_R(s_1;\mu )=\sum _{i\in s_2}({\pi }_{2i}^{-1}-1)e_i(\mu ,s_1{)}^2$$

is unbiased for $D_R(s_1;\mu )$ conditional on $s_1$ , since

$$E_s\{{\hat{D}}_R(s_1;\mu )|s_1\}=E_s\{\sum _{i\in s_2}{\pi }_{2i}^{-1}{e}_i(\mu ,s_1{)}^2-A_2|s_1\}=E_s\{D_R(s_1;\mu )|s_1\}.$$

(11)

Meanwhile, the TSE of the SRB-predictor $\overline{\mu }(x,s)$ given by Equation (5) is

$$D(s;\overline{\mu })=\sum _{i\in R}{e}_i(\overline{\mu }{)}^2\mathrm{and}{e}_i(\overline{\mu })=\overline{\mu }(x_i,s)-y_i.$$

For any $i\in R$ with $x_i=x$ , the errors of $\overline{\mu }(x,s)$ and $\mu (x,s_1)$ are related by

$$e_i(\mu ,s_1)=\mu (x,s_1)-y_i=\{\mu (x,s_1)-\overline{\mu }(x,s)\}+e_i(\overline{\mu }).$$

Since $\overline{\mu }(x,s)$ and $e_i(\overline{\mu })$ are constant of $s_1~q(s_1|s)$ , we have

$$e_i(\overline{\mu }{)}^2=E_q\{e_i(\mu ,s_1{)}^2|s\}-E_q\{a_i(\mu ,s_1{)}^2|s\}$$

(12)

where $a_i(\mu ,s_1)=\mu (x,s_1)-\overline{\mu }(x,s)$ and $E_q\{a_i(\mu ,s_1{)}^2|s\}$ is the variance of $\mu (x,s_1)$ under the SRB operation. Design-unbiased estimation of the risk $D(s;\overline{\mu })$ is given by the result below, the proof of which is given in Appendix.

• Theorem 2. For any given $\mu (·)$ , an unbiased estimator of the risk $\tau (\overline{\mu })$ of the corresponding SRB-predictor $\overline{\mu }(x,s)$ , over $s~p(s)$ , is given by

$$\hat{D}(s;\overline{\mu })=E_q(\sum _{i\in s_2}({\pi }_{2i}^{-1}-1)\{e_i{(\mu ,s_1)}^2-a_i{(\mu ,s_1)}^2\}|s).$$

In practice, where exact SRB is infeasible numerically, one can use the MC-SRB predictor based on $T$ subsamples, which is given as

$$\{\begin{array}{c}\tilde{\mu }(x_i,s)=T^{-1}\sum _{t=1}^T\mu (x_i,s_1^{(t)})\mathrm{if}i\in R\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\mu }(x_i,s)=T_i^{-1}\sum _{t=1}^T\mathbb{I}(i\notin s_1^{(t)})\mu (x_i,s_1^{(t)})\mathrm{if}i\in s}\end{array}$$

where $s_1^{(t)}$ is the $t$ -th subsample, $T_i=\sum _{t=1}^T\mathbb{I}(i\notin s_1^{(t)})$ and $s_2^{(t)}=s\backslash s_1^{(t)}$ .

To estimate the risk, for any $i\in s_2^{(t)}$ , let $e_i(\mu ,s_1^{(t)})=\mu (x_i,s_1^{(t)})-y_i$ directly, and let $a_i(\mu ,s_1^{(t)})=\mu (x_i,s_1^{(t)})-\stackrel{·}{\mu }(x_i,s)$ be an out-of-bag approximation to $\mu (x_i,s_1^{(t)})-\overline{\mu }(x_i,s)$ , instead of $\mu (x_i,s_1^{(t)})-\tilde{\mu }(x_i,s)$ that would have been a residual-based alternative. The MC risk estimator is given by

$$\tilde{D}(s;\overline{\mu })=\frac{1}{T}\sum _{t=1}^T\sum _{i\in s_2^{(t)}}({\pi }_{2i}^{-1}-1)\{e_i(\mu ,s_1^{(t)}{)}^2-a_i(\mu ,s_1^{(t)}{)}^2\}.$$

(13)

3.2. Using an Ensemble of Predictors

By design-based predictive inference, there is no need to assume that a true model exists for $y_U$ , or that one is able to identify the true model under repeated sampling. It is then natural to combine an ensemble of different predictors (e.g., Dietterich 2000; Zhou 2012; Dong et al. 2020; Sagi and Rokach 2018) in addition to selecting a single model and the corresponding predictor. Ensemble SRB prediction by voting or averaging is developed below.

3.2.1. SRB-Selector

Consider selecting a single model by voting given an order- $K$ heterogeneous ensemble $\{{\mu }_1,...,{\mu }_K\}$ . Let $D(s;{\mu }_k)=\sum _{i\in R}\{{\mu }_k(x_i,s)-y_i{\}}^2$ . Denote by $\Omega =\bigcup _{k=1}^K{\Omega }_k$ the partition of the sample space such that, for any $s\in {\Omega }_k$ and $l\ne k$ , we have

$$D(s;{\mu }_k)<D(s;{\mu }_l)$$

where we discount the possibility of $D(s;{\mu }_k)=D(s;{\mu }_l)$ merely to simplify the exposition. To select a single predictor for $R$ based on a given sample $s$ , which minimises the risk Equation (10), one would vote for ${\mu }_k(x,s)$ iff $s\in {\Omega }_k$ . The optimal selector is thus the perfect classifier of $\mathbb{I}(s\in {\Omega }_k)$ .

In practice it is a common approach to apply cross-validation and majority-vote, where cross-validation is based on $s_1~q(s_1|s)$ and $s_2=s\backslash s_1$ . The expected selection result is given by the SRB-selector as follows. Given any $(s_1,s_2)$ by $q(s_1|s)$ and any $k=1,...,K$ , let

$${\delta }_k(s_1)=\{\begin{array}{c}1\mathrm{if}k=\arg \underset{l=1,...,K}{min}\sum _{i\in s_2}{\{{\mu }_l(x_i,s_1)-y_i\}}^2\hfill \\ \multicolumn{1}{c}{0\mathrm{otherwise}}\end{array}$$

indicate which predictor has the least TSE in $s_2$ . The SRB-selector

$${\overline{\delta }}_k(s)=\{\begin{array}{c}1\mathrm{if}k=\arg \underset{l=1,...,K}{max}{E}_q\{{\delta }_k(s_1)|s\}\hfill \\ \multicolumn{1}{c}{0\mathrm{otherwise}}\end{array}$$

(14)

is a classifier of $\mathbb{I}(s\in {\Omega }_k)$ , that is, the expected majority-vote over cross-validation. Given the selection by Equation (14), say, ${\mu }_k$ , one can reuse the same cross-validation samples $(s_1,s_2)$ to obtain the selected SRB-predictor ${\overline{\mu }}_k$ and its associated risk.

3.2.2. Mixed SRB-Predictor

Consider averaging given an order- $K$ ensemble $\{{\mu }_1,...,{\mu }_K\}$ , and let the mixed SRB-predictor be

$$\mu (x,s)=\sum _{k=1}^K{w}_k{\overline{\mu }}_k(x,s)$$

(15)

where $\sum _{k=1}^K{w}_k=1$ and $w_k>0$ for the mixing weights, $k=1,\dots ,K$ . We have

$$D(s;\mu )=\sum _{k\ne 1}{w}_k^2{D}_{\mathrm{kk}}+{(1-\sum _{l\ne 1}{w}_l)}^2{D}_{11}+\sum _{k=1}^K\sum _{l\ne k,1}{w}_k{w}_l{D}_{\mathrm{kl}}+2\sum _{k\ne 1}{w}_k(1-\sum _{l\ne 1}{w}_l)D_{1k}$$

now that $w_1=1-\sum _{k\ne 1}{w}_k$ , where $D_{\mathrm{kk}}=D(s;{\overline{\mu }}_k)$ and $D_{\mathrm{kl}}=D(s;{\overline{\mu }}_k,{\overline{\mu }}_l)$ is given by

$$D_{\mathrm{kl}}=\sum _{i\in R}{e}_i({\overline{\mu }}_1)e_i({\overline{\mu }}_2)=\sum _{i\in R}{E}_q\{e_i({\mu }_1,s_1)e_i({\mu }_2,s_1)|s\}-E_q\{a_i({\mu }_1,s_2)a_i({\mu }_2,s_2)|s\},$$

that is, similarly to Equation (12). An estimator of $D_{\mathrm{kl}}$ follows as a corollary of Theorem 2, as well as its MC implementation similarly to Equation (13).

The optimal mixing weights $w_k$ minimise $D(s;\mu )$ . The estimated ${\hat{w}}_k$ can be obtained via $\hat{D}(s;\mu )$ given ${\hat{D}}_{\mathrm{kl}}$ , for all $k,l=1,...,K$ . Substituting ${\hat{w}}_k$ in Equation (15) yields the mixed SRB-predictor. The associated risk Equation (10) can be estimated by $\hat{D}(s;\mu )$ .

Whilst the above approach aims at minimising the risk, it may experience instability if the ensemble is not sufficiently heterogeneous. A robust approach to mixed ensemble prediction should automatically aim at the same mixing weight of two component predictors that are equal to other.

For any $k=1,...,K$ , write $\hat{D}(s;{\overline{\mu }}_k)=E_q({\hat{\tau }}_k|s)$ , similarly to $\hat{D}(s;\overline{\mu })$ in Theorem 2. Regarding the risk of ${\overline{\mu }}_k(x,s)$ defined by Equation (10), we have

$$\tau ({\overline{\mu }}_k)=E_{s_1}\{\tau ({\overline{\mu }}_k|s_1)\}=E_{s_1}\{E_s({\hat{\tau }}_k|s_1)\}=E_p\{E_q({\hat{\tau }}_k|s)\}$$

where $\tau ({\overline{\mu }}_k|s_1)$ is its conditional risk given $s_1$ . Let the SRB operation yield

$$w_k=E_q({\delta }_k|s),{\delta }_k=\{\begin{array}{c}1\mathrm{if}k=\arg \underset{l=1,...,K}{min}{\hat{\tau }}_l\hfill \\ \multicolumn{1}{c}{0\mathrm{otherwise}}\end{array}.$$

(16)

The corresponding mixed SRB-predictor, Equation (15), is robust against ${\mu }_k\approx {\mu }_l$ for any $k\ne l$ . While the SRB-selector, Equation (14), is a binary classifier taking the majority-vote over all $(s_1,s_2)$ , the robust mixing weight, Equation (16), is a proportion over all $(s_1,s_2)$ .

3.3. Illustration

Simulations below provide a simple illustration of design-based individual-level predictive inference. For better appreciation of the design-based approach, we include also the risk estimator under an assumed IID error model.

We generate 200 sets of $y_U$ of population size $N=2000$ in an ad hoc manner. For each $y_U$ , half of them are generated by M1 below and half of them by M2, where $x_1~N(0,1)$ and $x_2~\mathrm{Poisson}(5)$ ,

$$\begin{array}{c}\multicolumn{1}{c}{(\mathrm{M}1)y=x_1+0.5x_2+ϵ,ϵ~\{\begin{array}{c}N(0,1)\mathrm{if}z=1\iff x_2<3\hfill \\ \multicolumn{1}{c}{N(-2,1)\mathrm{if}z=2\iff 3\le x_2<7}\\ \multicolumn{1}{c}{N(2,1)\mathrm{if}z=3\iff x_2\ge 7}\end{array}}\\ \multicolumn{1}{c}{(\mathrm{M}2)y=0.5+1.5x_1+x_2+ϵ,ϵ~z^2+N(0,0.25),z~N(0,1),}\end{array}$$

From each population we draw a sample of size $n=200$ either by SRSWOR or Poisson Sampling. For Poisson Sampling, we set ${\pi }_i^{-1}\propto 1+1/\exp (\alpha +0.5y_i)$ and $\sum _{i\in U}{\pi }_i=n$ , where $\alpha \in \{1,-0.1,-1\}$ leads to the coefficient of variation of ${\pi }_i$ over $U$ , denoted by ${\mathrm{cv}}_{\pi }$ , to be about 15%, 30%, and 45%, respectively. This illustrates a situation where sample selection may cause issues for uncertainty assessment by the IID model.

Let an order-3 model ensemble contain linear regression, random forest, and support vector machine. Let the feature vector be $x=(x_1,x_2)$ in all the cases. We use a 70-30 random split for subsampling of $(s_1,s_2)$ and let $T=50$ for relevant MC-SRB operations such as Equation (13). We obtain thus the SRB-predictor as described in Equation (5) corresponding to each model.

For each SRB predictor, we estimate its standardized risk, Equation (10), $\tau /|R|$ , as described before, where we have ${\pi }_{2i}\equiv n_2/(N-n_1)$ under SRSWOR given $n_1=|s_1|$ and $n_2=|s_2|$ , and we use ${\varphi }_{2i}$ given by Equation (8) instead of ${\pi }_{2i}$ under Poisson Sampling. Note that if $\hat{\tau }$ is unbiased for $\tau$ over repeated sampling from a given population, then it is also unbiased for $D/|R|$ over all the 200 populations.

The average of the 200 true $D/|R|$ for each SRB predictor will be referred to as average true MSE in the results below, which is given by

$$\mathrm{MS}{\mathrm{E}}_{\mathrm{true}}=\frac{1}{200}\sum _{b=1}^{200}\frac{1}{N-n}\sum _{i\in R^{(b)}}\{\tilde{\mu }(x_i,s^{(b)})-y_i{\}}^2$$

where $s^{(b)}$ denotes the $b$ -th simulated sample and $R^{(b)}$ the corresponding out-of-sample units. The proposed MSE estimator, Equation (13), will be compared to two MSE estimators that rely on the IID error model (e.g., James et al. 2013): the residual-based estimator and the cross-validation (CrV) estimator. The average of the latter two estimators over the 200 populations are given as

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{MS}{\mathrm{E}}_{\mathrm{resid}}=\frac{1}{200}\sum _{b=1}^{200}\frac{1}{n}\sum _{i\in s^{(b)}}\{\tilde{\mu }(x_i,s^{(b)})-y_i{\}}^2}\\ \multicolumn{1}{c}{\mathrm{MS}{\mathrm{E}}_{\mathrm{CrV}}=\frac{1}{200T}\sum _{b=1}^{200}\sum _{t=1}^T\frac{1}{n_2}\sum _{i\in s_2^{(b,t)}}\{\mu (x_i,s_1^{(b,t)})-y_i{\}}^2}\end{array}$$

where $s^{(b)}=s_1^{(b,t)}\cup s_2^{(b,t)}$ signifies the $t$ -th subsampling of the $b$ -th sample.

Table 2 displays average true MSE and its estimates across the simulation settings. Regardless the model, the proposed design-based risk estimator Equation (13) is unbiased under SRSWOR $p$ -design where ${\pi }_{2i}$ is known. Whereas, using ${\varphi }_{2i}$ instead of ${\pi }_{2i}$ under informative Poisson sampling, it remains essentially unbiased when $\left({cv}_{\pi }=15\%\right)$ or $30\%$ , but the approximation may be seen to have caused some underestimation as ${\mathrm{cv}}_{\pi }$ increases to 45%, where the severest underestimate is $(9.288-9.866)/9.866\times 100=-5.86\%$ .

Table 2.

MSE and Estimates Given Each Model, Averaged Over 200 Simulations.

MSE	SRSWOR			PS $\left({cv}_{\pi }=15\%\right)$
	LR	RF	SVM	LR	RF	SVM
Average, true	8.399	9.013	9.272	8.566	9.225	9.671
Design, proposed	8.409	9.073	9.326	8.416	9.182	9.615
Model, CrV	8.457	9.481	9.862	8.014	9.214	9.405
Model, residual	8.162	5.105	7.706	7.766	4.945	7.578
MSE	PS $\left({cv}_{\pi }=30\%\right)$			PS $\left({cv}_{\pi }=45\%\right)$
	LR	RF	SVM	LR	RF	SVM
Average, true	8.957	9.726	10.451	9.866	10.884	11.573
Design, proposed	8.711	9.559	10.196	9.288	10.364	10.974
Model, CrV	7.624	8.880	8.799	6.992	8.262	7.933
Model, residual	7.369	4.731	7.330	6.776	4.367	6.758

Note. PS = Poisson sampling; LR = linear regression; RF = random forest; SVM = support vector machine.

The CrV-based IID-model MSE estimator is also essentially unbiased under SRSWOR, because the out-of-bag squared errors in the test sample $s_2$ have the same mean as those in $R$ conditional on $s_1$ under the $\mathrm{pq}$ -design. This is reasonable since the IID error model would hold exactly under SRS with replacement. As mentioned in Subsection 2.2, Bates et al. (2024) explain why the CrV-based MSE estimator is not exactly unbiased for the full-sample once-trained predictor $\mu (x,s)$ , even when the IID model is correct. The CrV-based MSE estimator is instead applied to the SRB-predictor $\overline{\mu }(x,s)$ here.

The CrV-based MSE estimator can become severely biased though, if the IID model does not hold for the actual sample selection mechanism, as illustrated here for Poisson Sampling with increasing ${\mathrm{cv}}_{\pi }$ . Furthermore, residual-based MSE estimation should be avoided even under SRSWOR $p$ -design, since it generally leads to large biases for predictors derived from highly flexible machine learning models such as random forest and support vector machine.

Next, to illustrate inference for ensemble individual prediction, we obtain the SRB-predictor Equation (5) selected by Equation (14), and the two mixed SRB-predictors using the weights that are either optimal for Equation (15) or robust Equation (16). Given any ensemble MC-SRB predictor $\tilde{\mu }(x,s^{(b)})$ , for $b=1,...,200$ , we calculate its design-based risk estimator using Equation (13) and following the corresponding description in Subsection 3.2. Whereas the two IID model-based MSE estimators are calculated in the same way as given above.

Table 3 presents average true MSE and estimates for the SRB selector and mixed SRB predictors across the simulations. The average true MSE by the SRB-selector is similar to that of the linear model in Table 2, where the MSE is the smallest by this model than random forest or support vector machine. The two mixed SRB predictors achieve largely the same true MSE for individual prediction, in each simulation setting, illustrating the robustness of ensemble prediction approach even when it cannot improve on the best single model in the given ensemble. In terms of MSE estimation, the results in Table 3 are seen to be consistent with what we have observed for Table 2.

Table 3.

MSE and Estimates for Ensemble Individual Prediction, Averaged Over 200 Simulations; Predictor Selected by Majority-Vote, or Averaged by Optimal or Robust Mixing Weights.

MSE	SRSWOR			PS $\left({cv}_{\pi }=15\%\right)$
	Selected	Optimal	Robust	Selected	Optimal	Robust
Average, true	8.432	8.367	8.380	8.570	8.558	8.578
Design, proposed	8.395	8.260	8.284	8.412	8.343	8.372
Model, CrV	8.453	8.341	8.374	8.015	7.981	8.009
Model, residual	8.076	7.264	7.178	7.746	7.146	7.008
MSE	PS $\left({cv}_{\pi }=30\%\right)$			PS $\left({cv}_{\pi }=45\%\right)$
	Selected	Optimal	Robust	Selected	Optimal	Robust
Average, true	8.980	8.985	9.012	9.897	9.915	9.981
Design, proposed	8.700	8.657	8.691	9.281	9.257	9.316
Model, CrV	7.631	7.618	7.647	7.006	6.997	7.035
Model, residual	7.289	6.902	6.667	6.989	6.977	7.008

Note. PS = Poisson sampling.

4. Illustrative Application

Here we describe an illustrative application of design-based prediction estimation for the Spanish Structural Business Survey (SSBS). The SSBS provides information about the main structural and economic characteristics of businesses, such as employed personnel, turnover, purchases, personnel expenses, taxes, and investments. The target population consists of businesses classified in one of the following economic sectors: industrial sector, commercial sector, and services sector.

We take the year 2020 for reference, where the SSBS population contained 2,615,811 business units and the estimated total turnover was 1,785 billion euros, which is related to the total value of market sales of goods and services to third parties during the reference year. The SSBS estimation is traditionally based on the HT-estimator. One of motivations of this study, which is directly related to the SSBS, is the need to investigate whether it is possible to reduce the SSBS sample size, by developing and introducing more efficient estimation approaches.

4.1. $\mathrm{pq}$ -Design for SSBS

The SSBS sample contains both fully surveyed business units and other units that are mainly imputed. For our purpose here, we shall only consider the sample of fully surveyed units, which are selected using a stratified random sampling design. There is as usual an exhaustive (i.e., take-all) stratum for the largest businesses, which will be excluded from the application below, since sampling error does not exist there. In addition, any stratum with only one or two sample units will be removed, because some variance smoothing techniques would be needed for these strata in practice, which have no direct relevance to the theory of design-based predictive inference.

A total of 9,681 strata are retained in this way, which contain altogether 2,018,561 population units. As shown in the top plot of Figure 2, the stratum population size is relatively small for most strata but has a skewed distribution. The biggest stratum does have 54,770 units and there are 319 strata with more than one thousand units. The histogram of stratum sample size is given in the bottom plot of Figure 2. Around a half of all the strata have five or fewer sample units, whereas the number of strata with sample size greater than 25 is 339. The total sample size is 80,280.

Figure 2.

Histogram of stratum size in population (top) and sample (bottom).

To investigate the potentials of sample size reduction, we selected randomly and without replacement 45% of the original sample units in each stratum, subjected to a minimum of three sample units in each stratum. The resulting total sample size is 40,514, which is about half of the original sample size. Alternative model-assisted estimators (to be described below) will be compared based on this reduced sample, which would demonstrate both the proposed inference approach and the potentials of sample reduction.

We notice that although the ad hoc reductions of stratum sample sizes above should not be taken as a proposal for the new SSBS sampling design, the estimation results based on this realized sample are more “tangible” than the alternative, whereby one first estimates the MSE of a given estimator based on the original sample and then speculates how this MSE might have changed had the total sample size been reduced by 50%. Moreover, insofar as our aim here is not the specifics of the new SSBS sampling design, the single realized sample is large enough to warrant the comparison of different estimation approaches, and there is no need to simulate the sample reduction above many times which would only have generated largely similar results.

We adopt within-stratum SRSWOR as the subsampling $q$ -design for any SRB-based estimator, in addition to the stratified random sampling $p$ -design above. A constant subsampling rate will be implemented in all the strata, where the stratum subsample sizes of $(n_1,n_2)$ are truncated to integers by the floor function, subjected to the constraint that $n_2\ge 2$ in each stratum, such that variance estimation is feasible. Although other subsampling $q$ -designs and stratum-specific subsampling rates can be explored, a systematic investigation in this respect is beyond the scope of this paper in any case.

4.2. Models and Estimators

We use turnover as the target variable for illustration. The HT-estimators based on the original SSBS sample (with 80,280 units) and the reduced sample (with 40,514 units) provide the baselines for comparison. Four additional estimators will be applied to the reduced sample, which arise from the $2\times 2$ combination of models (linear, tree) and estimators (prediction, unbiased) as described below.

The model is either linear or tree regression, defined globally regardless the design strata. The linear model uses four features:

• • the “administrative” turnover from the corporate incoming tax if available, or imputed stratum-mean (of available administrative turnover) if missing;

• • a binary indicator for whether the administrative turnover is missing;

• • the operating income according to the tax administration if available, or imputed stratum-mean (of available operating income) if missing;

• • a binary indicator for whether the operating income is missing.

The linear regression coefficients are estimated by weighted ridge regression, given the sampling design weights and a small tuning parameter $\lambda =0.01$ for the regularization penalty.

The tree regression model uses the following four features, where the missing values are not imputed but left as-is to the software package:

• • the administrative turnover,

• • the operating income,

• • the first digit of the National Classification of Economic Activities,

• • the number of employees according to the Business Register.

The tree model is built using the ready-made R package h2o for random forest, where one feature is chosen for each split (mtries = 1), the maximum tree depth is 20 (max_depth = 20) and the minimum number of observations per leaf is 5 (min_rows = 5). The observations are again weighted by the sampling weights.

Given either model, we consider the SRB prediction estimator Equation (4) and the design-unbiased model-assisted SRB-estimator (Sanguiao-Sande and Zhang 2021), where the latter can be given as ${\hat{Y}}_M=E_q({\hat{Y}}_{1M}|s)$ and

$${\hat{Y}}_{1M}=\sum _{i\in s_1}{y}_i+\sum _{i\in U\backslash s_1}\mu (x_i,s_1)+\sum _{i\in s_2}({\pi }_{2i}^{-1}-1)\{\mu (x_i,s_1)-y_i\}$$

The SRB operation uses a constant 80-20 sample-split for the linear model and a constant 50-50 sample-split for the tree model. Notice that in the case of tree model, the SRB-predictor Equation (5) is a random forest by construction since it is then the average of $T$ randomly constructed trees.

It is worth pointing out that, unlike the prediction estimator Equation (4) that applies $\overline{\mu }(x,s)$ directly to all the out-of-sample units, the SRB-estimator ${\hat{Y}}_M$ corrects the bias of each $\mu (x,s_1)$ using the observed errors $\{\mu (x_i,s_1)-y_i:i\in s_2\}$ via ${\pi }_{2i}$ ; see Sanguiao-Sande and Zhang (2021) for the details. Although ${\hat{Y}}_M$ is thus exactly design-unbiased, it may have a larger MSE than the prediction estimator Equation (4) that uses the same model. Conversely, a prediction estimator would become less attractive if its bias is ‘intolerable’, even though it may have a much smaller MSE than an unbiased estimator that uses the same model.

4.3. Results

Table 4 summarizes the results for the estimators described above. Apart from each estimate $\hat{Y}$ , it is also given the estimated bias (zero for an unbiased estimator), the estimated MSE, the relative error (RErr) given as $\sqrt{\mathrm{MSE}}/\hat{Y}$ , and the Monte Carlo (MC) error of the MSE estimator.

Table 4.

Estimation Results (in Billion Euros), $T=10,000$ Sample-Splits, Based on Same Reduced Sample Size Unless Indicated Otherwise.

Estimator, model	$\hat{Y}$	Bias	MSE	RErr	MC error
HT-estimator (full sample size)	258	0	94	0.04	—
HT-estimator	252	0	151	0.05	—
SRB-prediction estimator, linear	227	−2	50	0.03	3
SRB-prediction estimator, tree	238	4	27	0.02	5
SRB-estimator, linear	229	0	122	0.05	1
SRB-estimator, tree	234	0	107	0.04	2

Exact Rao-Blackwellisation for MSE estimation is simply beyond reach in this case, where we have more than six thousand strata with sample size three in the reduced sample. If we leave out two units in each stratum under subsampling, then there are more than $3^{6000}$ distinct samples under the $q$ -design just for these strata, because the models are not built separately within each stratum, not to mention the other strata with more sample units.

The MC error is the bootstrap estimated standard deviation of the MC-MSE estimator, which is due to the loss of MSE estimation efficiency by Monte Carlo compared to exact Rao-Blackwellisation (with zero MC error). It can be seen that the relative MC error still does not vanish despite the large number of sample-splits $T=10,000$ , for example, it is $3/50$ for the linear-model SRB-prediction estimator, and not surprisingly, the relative MC error can increase as the MSE reduces, such as $5/27$ for the tree-model SRB-prediction estimator. Nevertheless, the MSE estimation results here are reliable enough for us to distinguish between the different estimators.

First, we notice that the two design-unbiased SRB-estimators have only led to moderate efficiency gains over the HT-estimator. The main reason is likely to be the large number of strata with very few sample units. Basically, in case the SRB-estimator uses $\mu (x)$ that is given in advance, it would become a stratified difference estimator, which corrects for the design-based bias of $\mu (x)$ for each stratum population total by using only the within-stratum sample units. The situation is largely similar for the SRB-estimator using $\mu (x,s)$ estimated based on the whole sample, which may have a small sampling variance itself.

In comparison, the SRB-prediction estimators using the same models can be much more efficient precisely because they do not apply bias correction, that is, they are no longer stratified estimators as the SRB-estimators are, such that they can take full advantage of the reduced variance if the prediction biases are small. The MSE of the linear-model SRB-prediction estimator is only about one third of the sampling variance of the HT-estimator, given the same sample size, whereas the tree model further reduces the MSE by about 50% compared to the linear model. Meanwhile, while the bias of the linear-model prediction estimator is relatively small compared to its root MSE, this is no longer the case for the tree-model prediction estimator, that is, 4 against $\sqrt{27}$ .

This serves to remind one that it is often possible to reduce the MSE at some cost of increasing bias, such as when adopting either model-assisted or model-based estimators traditionally. The theory of design-based predictive inference allows one to estimate both the bias and MSE of a large class of prediction estimators Equation (1). This increases the scope of choice in practice, in order to achieve a sensible trade-off between bias and variance.

Finally, the illustrative results above suffice to demonstrate the potentials of sample size reduction for the SSBS. An appropriate scheme of stratum sample size reduction requires a more systematic investigation though. In particular, the accompanying estimator can be chosen from the broad class of prediction estimators Equation (1), assisted by any model or algorithm and the various features available from the administrative source. But a detailed analysis is needed to take into account the level of dissemination (instead of just an overall total here) and the tolerable trade-off between bias and root MSE.

5. Final Remarks

We have developed a theory of design-based predictive inference from finite-population probability sampling. For population total estimation, one would be interested in the total of the out-of-sample prediction errors, whereas the risk of individual-level prediction depends on the out-of-sample squared prediction errors. The SRB approach provides a unified treatment of both.

Adopting design-based predictive inference for official statistics allows one to circumvent the design versus model controversy. In addition to producing population-level estimates, it provides a theoretical basis for producing statistical registers or census-like data for descriptive official statistics. The theory we propose allows for any assisting ML models or algorithms, which can be more efficient than calibration estimation using only auxiliary totals or the parametric assisting models commonly used in survey sampling.

There are a number of issues worth further investigation, of which we only mention a few here. First, survey nonresponse is unavoidable in practice. Lee et al. (2023) apply a related SRB ensemble learning approach to missing data imputation. It would be helpful to develop a unified SRB approach, which can incorporate survey response under an extended quasi-randomisation framework. Next, other choices of risk than the total squared prediction errors may be considered for individual prediction, or interval estimation may be developed for population total inference wherever the design-based bias of the prediction estimator is deemed non-negligible. Finally, it is worth studying how better to balance between the risk of individual prediction and the MSE of population total estimation associated with the prediction estimator Equation (1).