A method to find an efficient and robust sampling strategy under model uncertainty

1 Introduction

We consider the problem of choosing strategy, in particular the design, for the estimation of the total of a study variable in a finite population when a set of $J$ auxiliary variables is available in a list sampling frame. We focus on the estimation of the total.

The decision about sampling strategy involves parameters which are unknown at the stage when the decision needs to be taken. After data collection the parameters can be estimated, although sometimes only under some assumptions. In practice, we often use data from previous waves of a repeated survey, frame variables or data from another survey that is similar to the one that under a planning stage. There is a risk that the available data do not give reliable information about relevant parameters. The method presented here involves a risk measure, which takes into account the possibility of being misled by inaccurate or incorrect beliefs about the values of the needed parameters. The risk measure is derived for the difference and the generalized regression estimators. Other than that, the measure is general. This measure and the discussion of its practical use are the main result of this paper.

One aim when selecting and devising the sampling strategy is efficiency in terms of small mean-squared error. The definition of “efficiency” is not unique, however, as it depends on the inference approach. Under the design-based approach, Godambe (1955), Lanke (1973) and Cassel et al. (1977) show that there is no uniformly best linear estimator, in the sense of being best for all populations. There is no best design either. Therefore, a traditional approach for defining the strategy has been to assume that the finite population is a realization of some superpopulation model. The strategy is then defined in such a way that it minimizes the model expected value of the design mean-squared error, a parameter called anticipated mean-squared error. The adjective “anticipated” was first introduced by Isaki and Fuller (1982) to emphasize the fact that this is a conceptual mean-squared error which is calculated in advance to sampling, based only on information available prior to sampling.

Assuming that a superpopulation model holds and its parameters are known, several authors have shown that the optimal strategy should make use of a probability proportional-to-size sampling design (e.g. Hájek, 1959; Cassel et al., 1976; Nedyalkova and Tillé, 2008). In practice, however, there is not even a consensus about the existence of a generating model, let alone what model to rely on. And even if there is a model, its parameters are unknown. There is evidence, rather empirical, that probability proportional-to-size sampling is not robust towards model misspecifications (e.g. Holmberg and Swensson, 2001). A second result of this paper is to provide some theoretical evidence of this fact.

Many articles discuss robustness in the survey sampling field. Beaumont et al. (2013), for instance, propose a robust estimator that downweights influential observations; Royall and Herson (1973) consider robustness under polynomial models; Bramati (2012) and Zhai and Wiens (2015) propose robust stratification methods. We provide theoretical evidence of lack of robustness of proportional-to-size sampling and propose a method for assisting in the decision about the sampling design.

The contents of the paper are arranged as follows. The optimal strategy under the superpopulation model is defined in section 2. The lack of robustness of this strategy when the model is misspecified is studied in section 3. The method for assisting on the choice of the sampling design is presented in section 4. In section 5, the risk measure introduced in the previous section is extended to be used together with the GREG estimator. Section 6 presents numerical illustrations of the results in the paper. First, we illustrate the lack of robustness of probability proportional-to-size sampling and the flexibility of the GREG estimator with a small simulation study. Second, we illustrate the implementation of the risk measure with real survey data. Finally, section 7 presents some conclusions.

2 Optimal strategy under the superpopulation model

Let $U$ be a finite population of size $N$ with elements labeled ${1, 2, \dots, k, \dots, N}$ . Let $x_{k} = (x_{1 k}, x_{2 k}, \dots, x_{J k})$

be a known vector of values of

J

auxiliary variables and

y_{k}

the unknown value of a study variable associated to unit

k \in U

. We are interested in the estimation of the total of

y

t_{y} = \sum_{U} y_{k}

Let $Ω$ be the power set of $U$ . A sample is any subset $s \in Ω$ and a sampling design

is a probability distribution on

Ω

, denoted by

P (S = s)

or simply

p (s)

. Let

π_{k} = \sum_{s ∋ k} p (s)

be the inclusion probability of

k

and

π_{k l} = \sum_{s \supset {k, l}} p (s)

the joint inclusion probability of

k

and

l

. A probability sampling design is a sampling design such that

π_{k} > 0

for all

k \in U

An estimator is a real valued function of the sample, ${^t}_{y} = {^t}_{y} (S)$ . By strategy we refer to the couple sampling design and estimator, $(p (\cdot), {^t}_{y})$ .

We consider only probability sampling designs with fixed sample size. As a convenient stepping stone we begin by considering unbiased linear estimators of the form

{^t}_{y} = (\sum U z_{k} - \sum s \frac{z_{k}}{π_{k}}) + \sum s \frac{y_{k}}{π_{k}} = \sum U z_{k} + \sum s \frac{e_{k}}{π_{k}}

(2.1)

with $z_{k}$ arbitrary known constants and $e_{k} = y_{k} - z_{k}$ . This estimator is called the difference estimator. The estimator defined in this way is said to be calibrated on $z$ as it satisfies ${^t}_{z} = \sum_{U} z_{k}$ . Note that if $z_{k} = 0$ for all $k \in U$ the estimator reduces to ${^t}_{y} = \sum_{s} y_{k} / π_{k}$ , that is, the Horvitz-Thompson estimator (Horvitz and Thompson, 1952). In later sections we focus on the generalized regression estimator (GREG).

The design MSE of the difference estimator is

{MSE}_{p} ({^t}_{y}) = {MSE}_{p} (\sum s \frac{e_{k}}{π_{k}}) = \sum U \sum U (π_{k l} - π_{k} π_{l}) \frac{e_{k}}{π_{k}} \frac{e_{l}}{π_{l}} .

(2.2)

As mentioned in the introduction, due to the non-existence of an optimal strategy under the design-based approach, often a superpopulation model, $ξ_{0}$ , is proposed and we search for an optimal strategy with respect to the anticipated mean-squared error,

{MSE}_{ξ_{0} p} ({^t}_{y}) = E_{ξ_{0}} {MSE}_{% p} ({^t}_{y}) = E_{ξ_{0}} E_{p} (({^t}_{y} - t_{y})^{2})

(2.3)

We may assume that the $y$ -values are realizations of the following model, denoted $ξ_{0}$ ,

Y_{k} = f (x_{k} | δ_{1}) + ϵ_{k} with E_{ξ_{0}} (ϵ_{k}) = 0, V_{ξ_{0}} (ϵ_{k}) = σ_{0}^{2} g (x_{k} | δ_{2})^{2} % and E_{ξ_{0}} (ϵ_{k} ϵ_{l}) = 0 \forall k \neq l

(2.4)

where $δ = (δ_{1}, δ_{2})$ is a vector of parameters, $f : R^{J} ⟶ R$ and $g : R^{J} ⟶ R^{+}$ . Following Rosén (2000a), the terms $f (x_{k} | δ_{1})$ and $g (x_{k} | δ_{2}) > 0$ will be called trend and spread, respectively. The term trend should not in general be understood in a temporal sense, rather it refers to the development of $y$ -values with $x$ .

Note that under $ξ_{0}$ , $e_{k}$ in the difference estimator (2.1

) is a random variable that represents the distance between the value of the study variable and

z_{k}

, i.e.

e_{k} = f (x_{k} | δ_{1}) + ϵ_{k} - z_{k}

. Therefore

E_{ξ_{0}} e_{k} = f (x_{k} | δ_{1}) - z_{k}

and

E_{ξ_{0}} e_{k}^{2} = (f (x_{k} | δ_{1}) - z_{k})^{2} + σ_{0}^{2} g (x_{k} | δ_{2})^{2}

. With some algebra, it can be seen from (2.2) and (2.3) that the anticipated MSE of the difference estimator becomes

{MSE}_{ξ_{0} p} ({^t}_{y}) = {MSE}_{p} (\sum s \frac{f (x_{k} | δ_{1}) - z_{k}}{π_{k}}) + σ_{0}^{2} (\sum U (\frac{1}{π_{k}} - 1) g (x_{k} | δ_{2})^{2})

(2.5)

Nedyalkova and Tillé (2008) derive the anticipated MSE in a more general case.

Tillé and Wilhelm (2017) give the anticipated MSE of the Horvitz-Thompson estimator. The second term in (2.5) is the Godambe-Joshi lower bound (e.g. Särndal et al. 1992, p. 453).

The anticipated MSE in (2.5) is the sum of two positive terms. It is easy to see that if

the estimator is calibrated on $z_{k} = f (x_{k} | δ_{1})$

the first term vanishes and the anticipated MSE equals the Godambe-Joshi lower bound

{MSE}_{ξ_{0} p} ({^t}_{y}) = σ_{0}^{2} (\sum U (\frac{1}{π_{k}} - 1) g (x_{k} | δ_{2})^{2})

(2.6)

Furthermore, after imposing the fixed sample size restriction $\sum_{U} π_{k} = n$ , if

the design is such that $π_{k} \propto g (x_{k} | δ_{2})$ , denoted $π$ ps( $g$ ),

the second term is minimized and we obtain

{MSE}_{ξ_{0} p}^{opt} ({^t}_{y}) = σ_{0}^{2} ⎛ ⎝ \frac{1}{n} {(\sum U g (x_{k} | δ_{2}))}_{k}^{2} - \sum U g (x_{k} | δ_{2})^{2} ⎞ ⎠ .

Conditions 1 and 2 suggest the specific roles of the design and the estimator in the sampling strategy. The estimator should “explain” the trend in the calibration sense of condition 1. The design should “explain” the spread. A strategy that satisfies conditions 1 and 2 simultaneously will be called optimal. In the same sense, any estimator and any design satisfying, respectively, condition 1 and 2, will be called optimal.

3 Robustness under a misspecified model

If the finite population is a realization of the superpopulation model (2.4), and if $f$ , $g$ and $δ$ were known, then an optimal strategy could be defined. In this section we study the robustness of this strategy when the model is misspecified.

We begin by defining how “misspecification”shall be understood in this paper. The working model $ξ_{0}$ reflects the beliefs the statistician has about the relation between the auxiliary variables $x$ and the study variable $y$ at the design stage. We shall assume that a true, unknown model $ξ$ exists. Any deviation of $ξ_{0}$ with respect to $ξ$ is a misspecification of the model. In order to keep the analysis tractable, we limit ourselves to the situation when the working model is of the form (2.4) and the true model, $ξ$ , is

Y_{k} = f (x_{k} | β_{1}) + ϵ_{k} with E_{ξ} (ϵ_{k}) = 0, V_{ξ} (ϵ_{k}) = σ^{2} g (x_{k} | β_{2})^{2} and E_{ξ} (ϵ_{k} ϵ_{l}) = 0 \forall k \neq l

(3.1)

where $β = (β_{1}, β_{2})$ is a vector of parameters, $f$ and $g$ as in (2.4) and $β \neq δ$ .

Result 1.

If $ξ_{0}$ is assumed when $ξ$ is the true superpopulation model, the model expected value of the design MSE in (2.2), under the difference estimator satisfying condition 1 above, becomes

{MSE}_{ξ p} ({^t}_{y}) = {MSE}_{p} (\sum s \frac{f (x_{k} | β_{1}) - f (x_{k} | δ_{1})}{π_{k}}) + σ^{2} (\sum U \frac{g (x_{k} | β_{2})^{2}}{π_{k}} - \sum U g (x_{k} | β_{2})^{2})

(3.2)

The result is proven by noting that $f (x_{k} | δ_{1})$ takes the role of $z_{k}$ in (2.5) and by taking into account that $e_{k} = f (x_{k} | β_{1}) - f (x_{k} | δ_{1}) + ϵ_{k}$ , therefore $E_{ξ} e_{k} = f (x_{k} | β_{1}) - f (x_{k} | δ_{1})$ and $E_{ξ} e_{k}^{2} = (f (x_{k} | β_{1}) - f (x_{k} | δ_{1}))^{2} + σ^{2} g (x_{k} | β_{2})^{2}$ . As the model is misspecified, we have deliberately avoided the use of the adjective “anticipated” in result 1.

Using result 1, it can be seen that for a design that satisfies condition 2 we obtain

{MSE}_{ξ, π ps} ({^t}_{y}) = {(\frac{\sum_{U} g (x_{k} | δ_{2})}{n})}^{2} {MSE}_{π ps} (\sum s \frac{f (x_{k} | β_{1}) - f (x_{k} | δ_{1})}{g (x_{k} | δ_{2})}) + σ^{2} (\frac{\sum_{U} g (x_{k} | δ_{2})}{n} \sum U \frac{g (x_{k} | β_{2})^{2}}{g (x_{k} | δ_{2})} - \sum U g (x_{k} | β_{2})^{2}) .

(3.3)

It is now possible to see that, even under a mild misspecification as the one considered here, the so-called optimal strategy is not necessarily optimal anymore, as its MSE (3.3) can be greater than the MSE obtained under other designs (3.2).

4 Guiding the choice of sampling design with the help of a risk measure

We have seen in section 3 that even a simple misspecification of the working model might result in the so-called optimal strategy not being optimal. It is therefore risky to accept a given model as correct without any type of assessment. While most of the information needed for an “objective” evaluation of the model is not available at the design stage, it is possible to reach some degree of confidence about the parameters in the working model that allows for comparing a set of designs and make the decision about which one to implement. In this section we propose a method to assist in the choice of the sampling design.

The expected MSE (3.2) in result 1 can be viewed as a function of $β$ and $σ^{2}$ , as everything else is available at the design stage. To begin with, let us assume that $σ^{2}$ is also known. Then we can write

L_{p} (β) = {MSE}_{ξ p} (β | x, δ, σ) = {MSE}_{p} (\sum s \frac{f (x_{k} | β_{1}) - f (x_{k} | δ_{1})}{π_{k}}) + σ^{2} (\sum U \frac{g (x_{k} | β_{2})^{2}}{π_{k}} - \sum U g (x_{k} | β_{2})^{2})

For any design, $p (\cdot)$ , this function can be evaluated at any $β$ and it indicates the loss incurred by assuming that $δ$ is the right parameter when it is, in fact, $β$ . We can assume a prior distribution on $β$ , $h (β)$ , and calculate the risk under $h$ ,

R (p) = E_{h} ({MSE}_{ξ p} (β | x, δ, σ)) = \int_{Θ} h (β) \cdot {MSE}_{ξ p} (β | x, δ, σ) d β,

(4.1)

where $Θ$ is the sample space of $β$ . The design that yields the smallest risk shall be chosen.

In practice, $σ^{2}$ is unknown. We propose three ways for dealing with it. The first one is to redefine $β$ as $(β, σ^{2})$ and calculate the risk as above. The second one is to provide some “guess” about it. The third one is to take into account that (proof in the appendix)

σ^{2} \approx \frac{S_{f, f}}{^{2}} ⎛ ⎝ \frac{1}{R_{f, y}^{2}} - 1 ⎞ ⎠

(4.2)

where $S_{f, f} = \sum_{U} (f (x_{k} | β_{1}) - ¯ f)^{2} / N$ , $¯ f = \sum_{U} f (x_{k} | β_{1}) / N$ , $^{2} = \sum_{U} g (x_{k} | β_{2})^{2} / N$ and $R_{f, y}$ is the correlation between $f (x | β_{1})$ and $y$ . (In example 3 below, we give a more convenient expression in a special case.) Although $R_{f, y}$ is unknown, for repeated surveys we do have some previous knowledge about it. In other cases it is often possible to have some reasonable “guess” about it.

It remains to comment on the choice of the prior distribution $h (β)$ . The choice of the distribution and its parameters is subjective and defined by the statistician. Nevertheless, it should reflect the available knowledge about the model parameter $β$ . In particular, $h (β)$ should be centered around $β = δ$

. Its variance should reflect how confident we are about the working model. Note that a full confidence on the working model would be a density with all its mass at

β = δ

, in which case the risk (4.1) would be minimized by the

π

ps design given by condition 2 in section 2.

It might be argued that by introducing $h (β)$ an additional source of subjectivity has been added to the choice of the sampling design. The prior may add a certain Bayesian flavor to the process, but note that $h (β)$ is only needed for choosing the design. Hence, the inference is still design-based. Furthermore, relying on an assumed model is also subjective in choice of assumption and it does involve a risk. The risk measure in (4.1) allows for quantification of that risk.

5 The risk measure under the Generalized Regression Estimator

The difference estimator (2.1) requires that $δ_{1}$ is fully specified in order to calculate $f (x_{k} | δ_{1})$ , which is undesirable from a practical standpoint. The generalized regression –GREG– estimator is an alternative that allows for the estimation of all or some of the components of $δ_{1}$ at the cost of introducing a small bias. In this section we adapt the material in sections 2 to 4 to strategies using the GREG estimator.

We define the generalized regression estimator in a more general way than in Särndal et al. (1992) as follows. Let $a_{k}$ ( $k = 1, \dots, N$ ) a weight defined by the statistician and $δ_{1} = (δ_{1}^{*}, δ_{1}^{* *})$ where $δ_{1}^{*}$ are fixed and $δ_{1}^{* *}$ are to be estimated. Let also

{^δ}_{1 s}^{* *} = {argmin}_{δ_{1}^{* *}} \sum s \frac{(y_{k} - f (x_{k} | δ_{1}))^{2}}{a_{k} π_{k}}

and ${^δ}_{1 s} = (δ_{1}^{*}, {^δ}_{1 s}^{* *})$ . The GREG estimator is

{^t}_{g r e g} = (\sum U f (x_{k} | {^δ}_{1 s}) - \sum s \frac{f (x_{k} | {^δ}_{1 s})}{π_{k}}) + \sum s \frac{y_{k}}{π_{k}}

(5.1)

An approximation to the design MSE of the GREG estimator is of the form (2.2) with $e_{k} = y_{k} - f (x_{k} | {^δ}_{1 U})$ where ${^δ}_{1 U} = (δ_{1}^{*}, {^δ}_{1 U}^{* *})$ and

{^δ}_{1 U}^{* *} = {argmin}_{δ_{1}^{* *}} \sum U \frac{(y_{k} - f (x_{k} | δ_{1}))^{2}}{a_{k}}

Example 1.

Let us consider the case where $f (x_{k} | δ_{1}) = δ_{1, 1} x_{1 k}^{δ_{1, J + 1}} + δ_{1, 2} x_{2 k}^{δ_{1, J + 2}} + \dots + δ_{1, J} x_{J k}^{δ_{1, 2 J}}$ . Let $δ_{1}^{*} = (δ_{1, J + 1}, \dots, δ_{1, 2 J})$ , $δ_{1}^{* *} = (δ_{1, 1}, \dots, δ_{1, J})^{'}$ and $x_{k}^{δ} = (x_{1 k}^{δ_{1, J + 1}}, \dots, x_{J k}^{δ_{1, 2 J}})$ . In this case we obtain

{^δ}_{1 s}^{* *} = {(\sum s \frac{x_{k}^{δ^{'}} x_{k}^{δ}}{a_{k} π_{k}})}^{- 1} \sum s \frac{x_{k}^{δ^{'}} y_{k}}{a_{k} π_{k}} and {^δ}_{1 U}^{* *} = {(\sum U \frac{x_{k}^{δ^{'}} x_{k}^{δ}}{a_{k}})}^{- 1} \sum U \frac{x_{k}^{δ^{'}} y_{k}}{a_{k}} .

Letting the exponents $δ_{1}^{*} = (δ_{1, J + 1}, \dots, δ_{1, 2 J}) = (1, \dots, 1)$ , we obtain the classical expression of the GREG estimator found in Särndal et al. (1992).

Example 2.

The case with only one auxiliary variable, i.e. $f (x_{k} | δ_{1}) = δ_{10} + δ_{11} x_{k}^{δ_{12}}$ with $a_{k} = 1$ , $δ_{1}^{*} = δ_{12}$ and $δ_{1}^{* *} = (δ_{10}, δ_{11})^{'}$ is known as the regression estimator. In this case we obtain the well known result that the design MSE can be approximated by expression (2.2) with $e_{k} = y_{k} - f (x_{k} | {^δ}_{1 U})$ where $f (x_{k} | {^δ}_{1 U}) = {^δ}_{10} + {^δ}_{11} x_{k}^{δ_{12}}$ and

{^δ}_{11} = \frac{N \sum_{U} x_{k}^{δ_{12}} y_{k} - \sum_{U} x_{k}^{δ_{12}} \sum_{U} y_{k}}{N \sum_{U} x_{k}^{2 δ_{12}} - {(\sum_{U} x_{k}^{δ_{12}})}_{U}^{2}} and {^δ}_{10} = \frac{1}{N} \sum U y_{k} - {^δ}_{11} \frac{1}{N} \sum U x_{k}^{δ_{12}} .

The misspecified model

Let us consider again the situation where the statistician uses the working model (2.4) but the true model is of the form (3.1) with $β_{1} = (β_{1}^{*}, β_{1}^{* *})$ , where $β_{1}^{*}$ is the counterpart of the fixed component $δ_{1}^{*}$ . The following result states a condition under which result 1 is valid for the GREG estimator.

Result 2.

If ${^δ}_{1 U}^{* *}$ converges in distribution to some $δ_{1}^{* *}$ then

{MSE}_{ξ p} ({^t}_{g r e g}) \to {MSE}_{p} (\sum s \frac{f (x_{k} | β_{1}) - f (x_{k} | δ_{1})}{π_{k}}) + σ^{2} (\sum U \frac{g (x_{k} | β_{2})^{2}}{π_{k}} - \sum U g (x_{k} | β_{2})^{2})

(5.2)

where $δ_{1} = (δ_{1}^{*}, δ_{1}^{* *})$ .

The result is proven by using the fact that if $X d \to C_{1}$ and $Y d \to C_{2}$ then $X \cdot Y d \to C_{1} \cdot C_{2}$ .

Example 3 (Continuation of Example 1).

Let the working model be as in example 1 and the true model be $f (x_{k} | β_{1}) = β_{1, 1} x_{1 k}^{β_{1, J + 1}} + β_{1, 2} x_{2 k}^{β_{1, J + 2}} + \dots + β_{1, J} x_{J k}^{β_{1, 2 J}}$ . Let also $β_{1}^{*} = (β_{1, J + 1}, \dots, β_{1, 2 J})$ , $β_{1}^{* *} = (β_{1, 1}, \dots, β_{1, J})^{'}$ and $x_{k}^{β} = (x_{1 k}^{β_{1, J + 1}}, \dots, x_{J k}^{β_{1, 2 J}})$ . In this case, ${^δ}_{1 U}^{* *} d \to A β_{1}^{* *}$ , where

A = {(\sum U \frac{x_{k}^{δ^{'}} x_{k}^{δ}}{a_{k}})}^{- 1} \sum U \frac{x_{k}^{δ^{'}} x_{k}^{β}}{a_{k}},

and (5.2) becomes

{MSE}_{ξ p} ({^t}_{g r e g}) \to {MSE}_{p} ⎛ ⎝ \sum s \frac{(x_{k}^{β} - x_{k}^{δ} A) β_{1}^{* *}}{π_{k}} ⎞ ⎠ + σ^{2} (\sum U \frac{g (x_{k} | β_{2})^{2}}{π_{k}} - \sum U g (x_{k} | β_{2})^{2}) .

(5.3)

Example 4 (Continuation of Example 2).

Let the working model be as in example 2 and the true model be $f (x_{k} | β_{1}) = β_{10} + β_{11} x_{k}^{β_{12}}$ with $β_{1}^{*} = β_{12}$ and $β_{1}^{* *} = (β_{10}, β_{11})^{'}$ . It can be shown that (5.2) becomes

{MSE}_{ξ p} ({^t}_{g r e g}) \to β_{11}^{2} {MSE}_{p} (\sum s \frac{v_{k}}{π_{k}}) + σ^{2} (\sum U \frac{g (x_{k} | β_{2})^{2}}{π_{k}} - \sum U g (x_{k} | β_{2})^{2})

(5.4)

with

v_{k} = (x_{k}^{β_{12}} -^{β_{12}}) - (x_{k}^{δ_{12}} -^{δ_{12}}) \frac{S_{β, δ}}{S_{δ, δ}},

(5.5)

and

	$^{β_{12}}$	$= \frac{1}{N} \sum U x_{k}^{β_{12}}$	$S_{β, δ}$	$= \frac{1}{N - 1} \sum U (x_{k}^{β_{12}} -^{β_{12}}) (x_{k}^{δ_{12}} -^{δ_{12}})$
	$^{δ_{12}}$	$= \frac{1}{N} \sum U x_{k}^{δ_{12}}$	$S_{δ, δ}$	$= \frac{1}{N - 1} \sum U (x_{k}^{δ_{12}} -^{δ_{12}})^{2}$

Note that (5.4) does not depend on $β_{10}$ .

If the condition for result 2 holds, i.e. if ${^δ}_{1 U}^{* *}$ converges in distribution to some $δ_{1}^{* *}$ , and this value can be calculated or approximated at the design stage for any $β$ , then the risk (4.1) can be computed.

Clearly, for the model that has been used in examples 1 and 3, the risk (4.1) can be obtained by means of expression (5.3). Furthermore, for the particular case developed in examples 2 and 4, where $f (x_{k} | β) = β_{10} + β_{11} x_{k}^{β_{12}}$ and $f (x_{k} | δ) = δ_{10} + δ_{11} x_{k}^{δ_{12}}$ , an alternative approximation of $σ^{2}$ is (proof in the appendix)

σ^{2} \approx β_{11}^{2} F_{0} with F_{0} = \frac{1}{^{2 β_{2}}} \frac{S_{1, β}^{2}}{S_{1, 1}} (\frac{1}{R_{x, y}^{2}} - \frac{1}{R_{1, β}^{2}})

(5.6)

where

^{2 β_{2}} = \frac{1}{N} \sum U x_{k}^{2 β_{2}} S_{1, β} = \frac{1}{N} \sum U (x_{k} - ¯ x) (x_{k}^{β_{12}} -^{β_{12}}) S_{1, 1} = \frac{1}{N} \sum U (x_{k} - ¯ x)^{2}

with $| R_{x, y} | \leq | R_{1, β} |$ and $R_{1, β}$ and $R_{x, y}$ are, respectively, the correlation coefficients between $x$ and $x^{β_{12}}$ and between $x$ and $y$ . The latter is unknown but often some decent guess about it is available.

The approximation of $σ^{2}$ in (5.6) is more convenient than (4.2) as now we have that (5.4) is approximated by

{MSE}_{ξ p} ({^t}_{g r e g}) \approx β_{11}^{2} [{MSE}_{p} (\sum s \frac{v_{k}}{π_{k}}) + F_{0} (\sum U \frac{g (x_{k} | β)^{2}}{π_{k}} - \sum U g (x_{k} | β)^{2})]

with $v_{k}$ given by (5.5). This expression depends neither on the intercept $β_{01}$ nor the parameter $σ$ , and the slope $β_{11}$ becomes a proportionality constant that can be ignored.

6 Numerical examples

In sections 2 and 3 we have established that the strategy that couples $π$ ps sampling with the difference estimator is optimal under a superpopulation model, but it is not robust to misspecifications of this model. In subsection 6.1 we present a small Monte Carlo simulation study carried out to illustrate these results by comparing the optimal strategy and three alternatives.

In sections 4 and 5 we introduced a measure that allows for quantifying the risk of implementing a sampling design, so allowing to guide the choice of design. In subsection 6.2 we illustrate the use of the risk measure with real survey data.

6.1 Simulation study under a misspecified model

We compare the efficiency and robustness of four strategies through a simulation study. The four strategies to be compared are $π$ ps together with the difference estimator (which is optimal when the model is correct), $π$ ps together with the GREG estimator (optimal design), stratified simple random sampling –STSI– together with the difference estimator (optimal estimator) and STSI together with the GREG estimator.

Our implementation of $π$ ps makes use of Pareto $π$ ps (Rosén, 1997). There is a host of other schemes for drawing $π$ ps samples. Nevertheless, Pareto $π$ ps is a convenient method with good properties, see for example Rosén (2000a).

Our implementation of STSI makes use of model-based stratification (Wright, 1983). We consider $H = 5$ strata with boundaries defined using the cum $\sqrt{f}$ -rule on $g (x_{k} | δ_{2})$ as in Särndal et al. 1992 (p. 463) and the sample is allocated using Neyman allocation, $n_{h} \propto N_{h} S_{g h}$ . Using the cum $\sqrt{f}$ -rule may be suboptimal (see Särndal et al. 1992, p. 464) but the efficiency of stratification by a continuous size variable is fairly insensitive to the exact choice of boundaries.

We consider only misspecification of the spread. The trend term is of the form $f (x_{k} | β_{1}) = β_{10} + β_{11} x_{k}^{β_{12}}$ with $β_{10} = 1000$ , $β_{11} = 1$ and $β_{12} = 0.75$ , $1$ and $1.25$ . The true spread is $g (x_{k} | β_{2}) = x_{k}^{β_{2}}$ with $β_{2} = 0.5$ , $0.75$ and $1$ . The working spread is $g (x_{k} | δ_{2}) = x_{k}^{δ_{2}}$ with $δ_{2} = 0.5$ , $0.75$ and $1$ .

We will use the difference estimator (2.1) calibrated on $f (x_{k} | β_{1})$ . Regarding the GREG estimator, we will fix $β_{12}$ , whereas the coefficients $β_{10}$ and $β_{11}$ will be estimated.

The simulation is set out as follows. The population size is $N = 5000$ . The $x$

-values are independent realizations from a gamma distribution with shape

α = 4 / 100

and scale

λ = 1200

plus one unit, whereas

y_{k}

is a realization from a gamma distribution with shape and scale

α_{k} = \frac{(β_{10} + β_{11} x_{k}^{β_{12}})^{2}}{σ_{0}^{2} x_{k}^{2 β_{2}}} and λ_{k} = \frac{σ_{0}^{2} x_{k}^{2 β_{2}}}{β_{10} + β_{11} x_{k}^{β_{12}}},

where $σ^{2}$ was set in such a way that the correlation between $x$ and $y$ is $ρ = 0.95$ . The design MSE of a sample of size $n = 500$ is then computed for each strategy. The process is iterated $B = 5000$ times.

Correct model
$β_{12}$	$β_{2}$	$δ_{2}$	$π$ ps–dif	$π$ ps–GREG	STSI–dif	STSI–GREG
0.75	0.50	0.50	$2.78 \cdot 10^{5}$	99.9	57.3	57.3
0.75	0.75	0.75	$4.82 \cdot 10^{4}$	99.6	77.9	77.9
0.75	1.00	1.00	$1.90 \cdot 10^{4}$	99.1	83.3	83.3
1.00	0.50	0.50	$7.64 \cdot 10^{6}$	99.9	57.3	57.3
1.00	0.75	0.75	$7.20 \cdot 10^{5}$	99.7	77.9	77.9
1.00	1.00	1.00	$2.14 \cdot 10^{5}$	99.1	83.2	83.3
1.25	0.50	0.50	$1.46 \cdot 10^{8}$	99.9	57.3	57.3
1.25	0.75	0.75	$7.85 \cdot 10^{6}$	99.7	78.0	78.0
1.25	1.00	1.00	$1.81 \cdot 10^{6}$	99.2	83.2	83.3
Misspecified model
$β_{12}$	$β_{2}$	$δ_{2}$	$π$ ps–dif	$π$ ps–GREG	STSI–dif	STSI–GREG
0.75	0.50	0.75	$3.98 \cdot 10^{5}$	99.9	98.9	98.9
0.75	0.75	1.00	$6.45 \cdot 10^{4}$	99.5	114.5	114.5
0.75	1.00	0.50	$4.73 \cdot 10^{4}$	100.1	133.9	134.0
1.00	0.50	1.00	$2.14 \cdot 10^{7}$	99.9	185.7	185.7
1.00	0.75	0.50	$1.03 \cdot 10^{6}$	100.1	93.1	93.1
1.00	1.00	0.75	$2.77 \cdot 10^{5}$	99.8	89.0	89.1
1.25	0.50	0.75	$2.09 \cdot 10^{8}$	99.9	98.9	98.9
1.25	0.75	1.00	$1.05 \cdot 10^{7}$	99.6	114.6	114.6
1.25	1.00	0.50	$4.50 \cdot 10^{6}$	100.3	134.0	134.2

Table 6.1: Efficiency of three strategies as a percentage of the expected MSE of

π

ps–dif under a misspecified model

Table 6.1 shows the results of the simulation study. The first three columns indicate the model parameters. The fourth column shows the (simulated) expected MSE of the strategy $π$ ps–dif, whereas the last three columns show the (simulated) efficiency of the strategies $π$ ps–GREG, STSI–dif and STSI–GREG compared to $π$ ps–dif (as a percentage), with efficiency defined as

eff = \frac{1}{B} B \sum r = 1 {eff}^{(r)} where {eff}^{(r)} = \frac{{MSE}_{ξ, π ps}^{(r)} ({^t}_{y})}{% {M S E}_{ξ, p}^{(r)} ({^t}_{y})},

in such a way that a value of 100 indicates that the strategy is as efficient as $π$ ps–dif and values smaller (larger) than one indicate that the strategy is less (more) efficient than $π$ ps–dif.

The upper part of Table 6.1 shows the case when the working model coincides with the true model. As expected, the strategy that couples $π$ ps with the difference estimator ( $π$ ps–dif) was always more efficient than the remaining strategies. Nevertheless, the loss in efficiency due to estimating some parameters through the GREG estimator is negligible. On the other hand, there is a remarkable loss in efficiency due to the use of STSI instead of $π$ ps. Finally, it is noted from (2.6) that as the anticipated MSE for all strategies does not depend on the trend $f$ but only on the spread $g$ , the efficiency remains constant under the same value of $δ_{2}$ , independently of the value of $β_{12}$ .

The lower part of Table 6.1 shows some comparisons under a misspecified model, in particular, a misspecified spread. It can be noted that even under this mild misspecification of the model, $π$ ps–dif is not necessarily the best strategy anymore as the strategies using STSI were more efficient in several cases. However, it is not evident when will STSI be more efficient than $π$ ps or vice versa. The risk measure introduced in section 4 can be used to guide the choice between designs. The results shown in this section agree with those shown by for example Holmberg and Swensson (2001).

6.2 Using the risk measure for choosing the design in a real survey

In this subsection we illustrate the implementation of the risk measure using data from a real survey. We want to estimate $t_{y} = \sum_{U} y_{k}$ where $U$ is the set of residential properties in Bogotá, Colombia (of size $N = 681276$ ) and $y_{k}$ is the value of the $k$ th property in 2017 in COP. $x_{k}$ , the built-up area of the $k$ th property in square meters, is known for every $k \in U$ . The auxiliary variable $x$

has mean 184, standard deviation 110 and skewness 2.57. The desired sample size is

n = 1000

We assume that a model of the type $ξ_{0}$ with $f (x_{k} | δ_{1}) = δ_{10} + δ_{11} x_{k}^{δ_{12}}$ and $g (x_{k} | δ_{2}) = x_{k}^{δ_{2}}$ adequately describes the association between $x$ and $y$ . We plan to use the GREG estimator for estimating $δ_{10}$ and $δ_{11}$ , i.e. $δ_{1}^{* *} = (δ_{10}, δ_{11})$ . We will use the risk (4.1) in order to assist the decision between $π$ ps or STSI. $H = 6$ strata are used and we take $h (β_{12}, β_{2})$

as a bivariate normal distribution with no correlation between

β_{12}

and

β_{2}

. We consider two cases with different degrees of confidence regarding the working model.

Case 1.

In this case no information about $δ_{12}$ , $δ_{2}$ or $R_{x, y}$ is available. Naive values of $δ_{12} = 1$ , $δ_{2} = 1$ and $R_{x, y} = 0.75$ are considered. In order to reflect the uncertainty, $h (β)$ should have a large variance, therefore we set

[\begin{matrix} β_{12} β_{2} \end{matrix}] \sim N ([\begin{matrix} 1.0 1.0 \end{matrix}], [\begin{matrix} {0.3295}^{2} & 0 0 & {0.3295}^{2} \end{matrix}]) .

Evaluation of (4.1) yields $R (π ps) = 6.89 \cdot 10^{15} β_{11}^{2}$ and $R (st) = 1.59 \cdot 10^{15} β_{11}^{2}$ , suggesting that a stratified design should be used.

The design MSE of both strategies is computed and we get, ${MSE}_{π ps} ({^t}_{y}) = 2.29 \cdot 10^{25}$ and ${MSE}_{st} ({^t}_{y}) = 1.36 \cdot 10^{25}$ . The strategy suggested by (4.1) was indeed the best choice.

Case 2.

Using a sample from 2010, prior values of $δ_{12} = 1.9$ , $δ_{2} = 2$ and $R_{x, y} = 0.7$ are proposed. As the uncertainty here is smaller than that in case 1, we set a smaller variance,

[\begin{matrix} β_{12} β_{2} \end{matrix}] \sim N ([\begin{matrix} 1.9 2.0 \end{matrix}], [\begin{matrix} {0.2471}^{2} & 0 0 & {0.2471}^{2} \end{matrix}])

Evaluation of (4.1) yields $R (π ps) = 7.08 \cdot 10^{22} β_{11}^{2}$ and $R (st) = 4.06 \cdot 10^{18} β_{11}^{2}$ , suggesting that a stratified design should be used.

The design MSE of both strategies is computed and we get ${MSE}_{π ps} ({^t}_{y}) = 1.85 \cdot 10^{28}$ and ${MSE}_{st} ({^t}_{y}) = 1.91 \cdot 10^{25}$ . Note that the use of (4.1) prevented us from using $π$ ps, whose MSE is almost one thousand times bigger than the one under stratified sampling!

7 Conclusions

The strategy that couples $π$ ps with the difference estimator is optimal when the parameters of the superpopulation model are known. Taking into account that these assumptions are seldom satisfied, it was shown in section 3 and illustrated in subsection 6.1 that this optimality breaks down even under small misspecifications of the model.

In section 4 we propose a method for choosing the sampling design, which is extended to its use with the GREG estimator in section 5. The method allows for taking the uncertainty about the model parameters into account by introducing a prior distribution on them. Although it could be argued that a source of subjectivity is added by introducing a prior distribution on the parameters, our view is that it is more subjective to choose the design without any type of assessment of the assumptions. Furthermore, inference is still design-based, as the prior is used only for choosing the design.

The method was illustrated with a real dataset, yielding satisfactory results. It should be noted that although the illustrations used stratified simple random sampling, the method in this article is valid for any sampling design.

Appendix. Proof of (4.2)

Proof.

The following expectations are required in the proof,

	$E_{ξ} Y_{k}$	$= E_{ξ} [f (x_{k} \| β_{1}) + ϵ_{k}] = f (x_{k} \| β_{1})$		(.1)
	$E_{ξ} Y_{k}^{2}$	$= E_{ξ} [{(f (x_{k} \| β_{1}) + ϵ_{k})}_{k}^{2}] = f (x_{k} \| β_{1})^{2} + σ^{2} g (x_{k} \| β_{2})^{2}$		(.2)

$E_{ξ} ¯ ¯¯ ¯ Y$ , $E_{ξ}^{2}$ and $E_{ξ} ¯ ¯¯¯¯¯ ¯ f Y$ are obtained using (.1) and (.2),

	$E_{ξ} ¯ ¯¯ ¯ Y$	$= E_{ξ} [\frac{1}{N} \sum U Y_{k}] = \frac{1}{N} \sum U E_{ξ} Y_{k} = \frac{1}{N} \sum U f (x_{k} \| β_{1}) \equiv ¯ ¯ ¯ f$		(.3)
	$E_{ξ}^{2}$	$= E_{ξ} [\frac{1}{N} \sum U Y_{k}^{2}] = \frac{1}{} N \sum U (f (x_{k} \| β_{1})^{2} + σ^{2} g (x_{k} \| β_{2})$

A method to find an efficient and robust sampling strategy under model uncertainty

Authors

Adaptive Sampling for Stochastic Risk-Averse Learning

Active Sentence Learning by Adversarial Uncertainty Sampling in Discrete Space

Age-optimal Sampling and Transmission Scheduling in Multi-Source Systems

Adaptive Sampling for Convex Regression

Optimizing the CVaR via Sampling

Sub-sampling and other considerations for efficient risk estimation in large portfolios

Myopic robust index tracking with Bregman divergence

1 Introduction

2 Optimal strategy under the superpopulation model

3 Robustness under a misspecified model

Result 1.

4 Guiding the choice of sampling design with the help of a risk measure

5 The risk measure under the Generalized Regression Estimator

Example 1.

Example 2.

The misspecified model

Result 2.

Example 3 (Continuation of Example 1).

Example 4 (Continuation of Example 2).

6 Numerical examples

6.1 Simulation study under a misspecified model

6.2 Using the risk measure for choosing the design in a real survey

Case 1.

Case 2.

7 Conclusions

Appendix. Proof of (4.2)

Proof.

A method to find an efficient and robust sampling strategy under model uncertainty

Authors

Adaptive Sampling for Stochastic Risk-Averse Learning

Active Sentence Learning by Adversarial Uncertainty Sampling in Discrete Space

Age-optimal Sampling and Transmission Scheduling in Multi-Source Systems

Adaptive Sampling for Convex Regression

Optimizing the CVaR via Sampling

Sub-sampling and other considerations for efficient risk estimation in large portfolios

Myopic robust index tracking with Bregman divergence

This week in AI

1 Introduction

2 Optimal strategy under the superpopulation model

3 Robustness under a misspecified model

Result 1.

4 Guiding the choice of sampling design with the help of a risk measure

5 The risk measure under the Generalized Regression Estimator

Example 1.

Example 2.

The misspecified model

Result 2.

Example 3 (Continuation of Example 1).

Example 4 (Continuation of Example 2).

6 Numerical examples

6.1 Simulation study under a misspecified model

6.2 Using the risk measure for choosing the design in a real survey

Case 1.

Case 2.

7 Conclusions

Appendix. Proof of (4.2)

Proof.