11 Balanced Sampling

The cube method proposes a general procedure that allows the selection of balanced random samples, with equal or unequal inclusion probabilities, in the sense that the Horvitz-Thompson estimates are equal, or nearly equal, to the population total of the balancing variables. Till’e (2006)

Balanced sampling has commonly been known as a nonprobability sampling technique, such as quota sampling, convenience sampling, or judgment sampling. This type of sampling suggests selecting samples for which the sample mean of an auxiliary information characteristic is identical to the population mean of that auxiliary information characteristic. Moreover, if this auxiliary information characteristic is highly correlated with the characteristic of interest, then balanced sampling is said to be optimal because it will accurately reproduce the total or mean of the characteristic of interest in the population.

Till’e (2006) states that the idea of selecting balanced samples began with Neyman (1934), who said that “the method of purposive selection consists of a) dividing the population of districts into second-order strata according to the values of \(x\) and \(y\), and b) randomly selecting from each stratum a fixed number of districts. The number of selections is determined by the condition that the weighted average of the characteristic of interest be maintained.” Later, Yates (1946) includes the following excerpt: “A random sample must be selected. Individuals will be included through the same random process; the first member will be compared with the first member of the original sample, the second individual with the second member of the original sample, and so on. A new member will be substituted if it improves the balance.”

Recently, partial solutions have been obtained for the random selection of balanced samples, through properly defined sampling designs, using methods proposed by recognized authors such as Ardilly (1991) and Deville (1992). On the other hand, authors such as Valliant et al. (2000) and Royal and Herson (1973) have considered the construction of estimators, framed under inference methods based only on population models, and their optimality from the model point of view without considering the sampling design. They conclude that a sampling design can be balanced even if it is not necessarily random or probabilistic.

In contrast, Deville and Tillé (2004) developed a general and rigorous procedure that allows the extraction of balanced probability samples and the subsequent estimation of the quantities of interest, framed under design-based inference methods. This procedure is known as the cube method and allows the selection of random samples over a set of auxiliary information characteristics, or balancing variables. It has the attractive property that the Horvitz-Thompson estimator reproduces the population total of the balancing variables. Later, Deville and Till’e (2005) adapted a variance approximation for the Horvitz-Thompson estimator in balanced sampling.

11.1 Notation

Because, under a balanced sampling design, the Horvitz-Thompson estimator for the totals of a set of auxiliary variables must be equal to their population total, the variance of the estimator of the population total of the characteristic of interest should decrease as the correlation with the auxiliary variables increases.

The objective is to estimate the population total of the characteristic of interest, \(t_y=\sum_{k\in U}y_k\). Thus, it is assumed that the values of the vectors \[ \mathbf{x}_k=\left(x_{k1},x_{k2},\ldots,x_{kQ}\right)' \]

taken for \(q\) balancing variables are known for all units in the population. Therefore, the vector of totals of the balancing variables \[ \mathbf{t}_{\mathbf{x}}=\sum_{k \in U}\mathbf{x}_k' \]

is also known and can be estimated using the Horvitz-Thompson estimator through the following expression: \[ \hat{\mathbf{t}}_{\mathbf{x},\pi}=\sum_{k \in U}\frac{\mathbf{x}_k}{\pi_k}I_k. \]

The objective is to construct a balanced sampling design, defined as follows.

Definition

A sampling design is balanced with respect to the auxiliary variables \(x_1, ..., x_Q\) if and only if it satisfies the balancing equations given by \[ \hat{\mathbf{t}}_{\mathbf{x},\pi}=\mathbf{t}_{\mathbf{x}} \]

for every sample \(s \in \mathcal{S}\) such that \(p(\mathbf{s})>0\) and for every \(q=1,...,Q\). In other words, \[ Var(\hat{\mathbf{t}}_{\mathbf{x},\pi})=\mathbf{0} \]

Note that \(Var(\hat{\mathbf{X}}_{\pi})\) is a variance-covariance matrix. In these terms, the balanced sampling design defines a support \(\mathcal{Q}\) given by \[ \mathcal{Q}=\left\{ \mathbf{I} \in \mathcal{S}| \sum_{k\in U}\frac{\mathbf{x}_k}{\pi_k}I_k=\mathbf{t}_{\mathbf{x}} \right\} \]

where \(\mathbf{I}=(I_1,\ldots,I_n)'\) is the inclusion vector for the elements in the sample and \(\mathcal{S}\) is the symmetric without-replacement support. To accept that a sampling design can be conditioned, the reader should be familiar with the definitions given in the first chapters of this text. In particular, note that, from Definition 2.1.5, the symmetric without-replacement support, which allows the definition of simple random sampling, among others, is also a conditioned support and is given by \[ \mathcal{S}_n=\left\{ \mathbf{s} \in \mathcal{S}| \sum_{k\in U}s_k=n \right\} \]

Similarly, the symmetric with-replacement fixed-size support, which allows the proper definition of the simple random with-replacement design, among others, is conditioned because \[ \mathcal{R}_n=\left\{ \mathbf{s} \in \mathcal{R} | \sum_{k\in U}s_k=n \right\} \]

11.1.1 Examples

Some examples are presented below. Although they are not useful in practice, they illustrate the objective of balanced sampling.

Example

Simple random sampling: this class of fixed-size sampling designs of size \(n\) is balanced on the variable \(x_k=\pi_k\), \(k\in U\). Indeed, \[ \sum_{k\in S}\frac{xk}{\pi_k}=\sum_{k\in S}1=n=\sum_{k\in U}\pi_k \]

Example

Stratification: suppose that in a population stratified into \(H\) strata (\(U_h\), \(h=1,...,H\), \(\#U_h=N_h\)), a simple random sample of size \(n_h\) is selected in each stratum. The design is balanced on the variables \[ \delta_{kh}= \begin{cases} 1 & \text{if unit *k* is in stratum h,}\\ 0 & \text{otherwise} \end{cases} \]

Since \[ \sum_{k\in S}\frac{\delta_{kh}}{\pi_k}=\sum_{k\in S}\delta_{kx}\frac{N_h}{n_h} =N_h=\sum_{k\in U}\delta_{kh} \]

In most practical problems, the balancing equations cannot be exactly satisfied. In other words, there is a rounding problem because the inverse of the inclusion probability is not an integer. For this reason, the objective is to construct a sampling design that satisfies the balancing equations exactly, if possible, or to find the best approximation if it is not. The rounding problem is negligible when the expected sample size is large.

11.2 The Cube Method

This method is composed of two phases, called the flight phase and the landing phase. In the first phase, for the restrictions to be satisfied exactly, the inclusion probabilities must be rounded to zero (0) or one (1). The landing phase consists of the appropriate handling of the rounding.

As we have seen, each vector \(\mathbf{s}\), in sampling without replacement, is a vertex of an N-cube, and the number of possible samples is the number of vertices of the N-cube. A sampling design with inclusion probability vector \(\boldsymbol{\pi}\) consists of assigning a probability to each vertex.

Geometrically, a sampling design consists of expressing the vector \(\boldsymbol{\pi}\) as a convex linear combination of the vertices of the N-cube. An algorithm can be viewed as a random path that reaches a vertex of the N-cube in such a way that the balancing equations are satisfied.

11.2.1 Flight Phase

This is a random walk that begins with an inclusion probability vector and remains in the intersection between the cube and the subspace restricted by the balancing equations. This random walk stops at a vertex of that intersection.

The objective of this phase is to randomly choose a vertex of \[K=\{[0,1]^N \cap Q\},\] where \(Q=\boldsymbol{\pi}+\ker\mathbf{A}\) and \(\mathbf{A}=(\check{\mathbf{x}}_1, ..., \check{\mathbf{x}}_N)\), in such a way that the balancing equations are reproduced satisfactorily. The landing phase is necessary only if the selected vector is not a vertex of the cube, and it consists of relaxing the restrictions as little as possible in order to select a sample, that is, a vertex of the cube.

Example

The flight phase transforms an inclusion probability vector into a vector of zeros and ones. \[ \boldsymbol{\pi}=\begin{pmatrix} 0.5 \\ 0.5 \\ 0.5 \\ 0.5 \end{pmatrix}\rightarrow\begin{pmatrix} 0.666 \\ 0.666 \\ 0.666 \\ 0 \end{pmatrix}\rightarrow\begin{pmatrix} 1 \\ 0.5 \\ 0.5 \\ 0 \end{pmatrix}\rightarrow\begin{pmatrix} 1 \\ 0 \\ 1 \\ 0 \end{pmatrix} \]

If there is a rounding problem, then some components cannot be converted to zero: \[ \boldsymbol{\pi}=\begin{pmatrix} 0.5 \\ 0.5 \\ 0.5 \\ 0.5 \\ 0.5 \end{pmatrix}\rightarrow\begin{pmatrix} 0.625 \\ 0 \\0.625 \\ 0.625 \\ 0.625\end{pmatrix}\rightarrow\begin{pmatrix} 0.5 \\ 0 \\ 0.5 \\ 1 \\ 0.5\end{pmatrix}\rightarrow\begin{pmatrix} 1 \\ 0 \\ 0.5 \\ 1\\ 0\end{pmatrix} \]

11.2.2 The Balanced Martingale

The general algorithm for carrying out the flight phase is developed using the following definition.

Definition

A discrete random process \(\boldsymbol{\pi}(t) = [\boldsymbol{\pi}_k(t)]\) in \(\mathbb{R}^N\), \(t= 0, 1, ...\), is called a balanced martingale for an inclusion probability vector \(\boldsymbol{\pi}\) and for the auxiliary variables \(x_1, ..., x_p\) if

\(\boldsymbol{\pi}(0) = \boldsymbol{\pi}\),
\(E [\boldsymbol{\pi}(t)|\boldsymbol{\pi}(t-1), ....,\boldsymbol{\pi}(0)] = \boldsymbol{\pi}(t-1)\), \(t = 1, 2, ...\)
\(\boldsymbol{\pi}(t) \in K = \{[0, 1]^N \cap (\boldsymbol{\pi} + \ker A)\)

11.2.3 Implementation of the Flight Phase

First, we initialize with \(\boldsymbol{\pi}(0) = \boldsymbol{\pi}\). Then, at stage \(t = 1, ...., T\),

We define a vector \(\mathbf{u}(t)=[u_k(t)] \neq 0\) such that

\(\mathbf{u}(t)\) is in the kernel of matrix A,
\(u_k(t) = 0\) if \(\pi_k(t)\) is an integer.

We calculate \(\lambda_1^*(t)\) and \(\lambda_2^*(t)\), the largest values such that \[ 0\leq \boldsymbol{\pi}(t) + \lambda_1^*(t)u(t)\leq1, \] \[ 0\leq \boldsymbol{\pi}(t) - \lambda_2^*(t)u(t)\leq1, \]
We choose \[ \boldsymbol{\pi}(t) = \begin{cases} \boldsymbol{\pi}(t-1)+ \lambda_1^*(t)\mathbf{u}(t) & \text{with probability } q_1(t)\\ \boldsymbol{\pi}(t-1)- \lambda_2^*(t)\mathbf{u}(t) & \text{with probability } q_2(t) \end{cases} \] where \[ q_1(t)=\lambda_2^*(t)/(\lambda_1^*(t)+\lambda_2^*(t)) \] and \[ q_2(t)=\lambda_1^*(t)/(\lambda_1^*(t)+\lambda_2^*(t)) \]

11.2.4 The Landing Phase

At the end of the first phase, the balanced martingale has reached a vertex of \(K\), which is not necessarily a vertex of \(C\). This vertex is denoted as \(\boldsymbol{\pi}^*=[\pi_k^*]=\boldsymbol{\pi}(T)\). Let \(q\) be the number of noninteger components in this vertex. If \(q=0\), the algorithm is complete. If \(q>0\), some restrictions cannot be satisfied rigorously.

Let \(U=\{k \in U|0 < \pi_k^* < 1\}\). The objective is to find a sampling design that yields a sample \(s^*\subset U^*\) such that \[\sum_{k\in S} a_k \approx \sum_{k\in U} a_k\pi_k^* = \sum_{k\in U}a_k\pi_k,\] with \(a_k=\check{\mathbf{x}}_k\) and \(s^*=s \cap U^*\).

This is solved by linear programming. Applying the simplex method, we have \[ \min_{p^*(\cdot)}\sum_{s^*\subset U^*}Cost(s)p^*(s), \]

subject to \[ \begin{aligned} \sum_{s^*\subset U}p(s^*)&=1\\ \sum_{s^*\ni k}p(s^*)&=\pi_k\\ 0\leq p(s^*)& \leq 1 \end{aligned} \]

where \(Cost(s)\) is the cost of the sample, which increases if the balancing equations given in the previous sections are not met. A sample is then selected with a sampling design \(p(\cdot)^*\). This program does not depend on the population size, but only on the number of balancing variables. If the number of auxiliary variables is very large, one auxiliary variable must be removed at the end of the flight phase. For this reason, it is important to order the balancing variables according to their correlation with the variables of interest.

11.2.5 Variance

Deville and Till’e (2005) proposed approximating the variance by assuming that the balanced sampling measure can be treated as conditional Poisson sampling. Thus, \[ Var(\hat{t}_{y,\pi})= Var(\hat{E}_{poisson})= \frac{N}{N - p}\sum_{k\in U} \frac{E^2_k}{\pi_k^2}\pi_k(1-\pi_k), \]

where \(E_k = y_k - \mathbf{x}'_k\mathbf{B}\).

Example

Note that the same function fulfilled by balanced sampling is fulfilled by the \(\pi\)PT sampling design because, by virtue of knowing a characteristic of interest and following Result 4.3.2, it is guaranteed that the estimator of the population total of the auxiliary information characteristic, \(\hat{t}_{x,\pi}\), reproduces the population total of the characteristic of interest, \(t_x\), with zero variance.

However, the \(\pi\)PT sampling design fulfills this function for one and only one auxiliary information characteristic. When the researcher has access to several auxiliary information characteristics simultaneously, \(\pi\)PT sampling is no longer useful. In this sense, it can be said, loosely speaking, that the balanced sampling design is a generalization of the \(\pi\)PT sampling design.

This example illustrates the computational procedure for achieving the final objective of selecting a balanced sample. The MU284 population (Särndal et al. 1992) will be used for this purpose. First suppose, without loss of generality, that the plan is initially to use a \(\pi\)PT sampling design; any other sampling design could also be used. Using the inclusionprobabilities function from the sampling package, the inclusion probabilities induced by this sampling design, with probability proportional to the auxiliary information characteristic P75, are obtained. Note that the sample size is 50 units.

library(sampling)
data(MU284)
pik <- inclusionprobabilities(MU284$P75, 50)
sum(pik)

[1] 50

Suppose that we want to obtain a balanced sample with respect to all auxiliary information characteristics given by P75, CS82, SS82, S82, ME84, and REV84. To do this, we include all observed population values of these balancing variables in a matrix. Next, we use the samplecube function to obtain a sample that is balanced with respect to all population totals of all balancing variables.

X <- as.matrix(MU284[, c("P75", "CS82", "SS82", "S82", "ME84", "REV84")])
s <- samplecube(X, pik, order = 1, comment = TRUE)


BEGINNING OF THE FLIGHT PHASE
The matrix of balanced variable has 6  variables and  284  units
The size of the inclusion probability vector is  284 
The sum of the inclusion probability vector is  50 
The inclusion probability vector has  281  non-integer elements
Step 1  


BEGINNING OF THE LANDING PHASE
At the end of the flight phase, there remain  6 non integer probabilities 
The sum of these probabilities is  2 
This sum is  integer
The linear program will consider  15  possible samples
The mean cost is  0.043 
The smallest cost is  0.0039 
The largest cost is  0.1 
The cost of the selected sample is 0.0039

QUALITY OF BALANCING
      TOTALS HorvitzThompson_estimators Relative_deviation
P75     8182                       8182 -0.000000000000056
CS82    2583                       2550 -1.295603899473805
SS82    6301                       6229 -1.137877152915336
S82    13500                      13290 -1.554492539420864
ME84  505226                     505421  0.038576506998984
REV84 874017                     866484 -0.861914494345551

Note that the output of this function is very informative. For this particular case, both the flight phase and the landing phase were needed. At the end of the flight phase, six individuals had probabilities that were neither zero nor one. Therefore, the cube method needs the landing phase to reach convergence. In addition to the comments for each phase of the cube method, this function also returns a table describing the quality of the procedure in terms of relative deviation. The reader should not overlook the quality of the balancing. It is simply extraordinary that such accuracy can be achieved with a sample of only 50 units.

11.3 Frame and Lucy

This chapter closes with the implementation of the cube method for selecting balanced samples. Suppose that the researcher knows the structural behavior of some characteristics of interest, namely Income and Number of employees. To select a balanced sample, first set the inclusion probabilities according to a simple random sampling design. As usual, insert the matrix of observations of the characteristics of interest into the samplecube function.

library(TeachingSampling)
data(BigLucy)
n <- 2000
N <- nrow(BigLucy)
pik <- rep(n / N, N)
X <- as.matrix(BigLucy[, c("Income", "Employees")])

s <- samplecube(X, pik, order = 1, comment = TRUE)


BEGINNING OF THE FLIGHT PHASE
The matrix of balanced variable has 2  variables and  85296  units
The size of the inclusion probability vector is  85296 
The sum of the inclusion probability vector is  2000 
The inclusion probability vector has  85296  non-integer elements
Step 1  


BEGINNING OF THE LANDING PHASE
At the end of the flight phase, there remain  2 non integer probabilities 
The sum of these probabilities is  0.24 
This sum is  non-integer
The linear program will consider  3  possible samples
The mean cost is  0.000021 
The smallest cost is  0.000001 
The largest cost is  0.000041 
The cost of the selected sample is 0.000022

QUALITY OF BALANCING
            TOTALS HorvitzThompson_estimators Relative_deviation
Income    36634733                   36657619              0.062
Employees  5391992                    5393863              0.035

For this particular case, the samplecube function, which implements the cube method, needed both the flight phase and the landing phase to reach convergence. The flight phase ended with 2 elements whose inclusion probabilities were neither zero nor one. However, after the landing phase, a balanced sample was selected. Once again, the quality of the balancing cannot go unnoticed.

After selecting the balanced sample, it is time to obtain the corresponding estimates. In general, it is possible to use the E.piPS function from the TeachingSampling package because the general framework of balanced sampling fits the characteristics that govern Horvitz-Thompson estimation.

sam <- (1:length(pik))[s == 1]
pik.s <- pik[sam]
sample_data <- BigLucy[sam, ]
target_variables <- data.frame(
  Income = sample_data$Income,
  Employees = sample_data$Employees,
  Taxes = sample_data$Taxes
)

E.piPS(target_variables, pik.s)

                   N     Income Employees     Taxes
Estimation     84742 36657619.3 5393863.2 1008092.1
Standard Error     0   490893.3   61832.5   31277.5
CVE                0        1.3       1.1       3.1
DEFF             NaN        1.0       1.0       1.0

The results returned by the function are optimal, in the sense that, in addition to obtaining estimates close to the population total for the characteristic of interest, they also maintain the population totals of the characteristics of interest in the sampling design.

11.3.1 Some Questions

Till’e (2006) answers some questions that arise directly with respect to how this new method works in practice:

Why not use calibration instead of balancing?\ Stratification is a special case of balanced sampling; post-stratification is a special case of calibration. In stratification and balancing, the weights are not random. This makes it a better strategy. Calibration has the advantage of requiring only knowledge of the population totals of the auxiliary variables, whereas balancing requires knowledge of the values of the auxiliary variables for all units in the population.
How accurate is the estimation approximation in balanced sampling?\ Deville and Tillé (2004) showed that under realistic regularity conditions in practical settings, \[ \left|\frac{\hat{t}_{x_q,\pi}-t_{xq}}{t_{xq}}\right|<O(p/N)\leq o_p(\sqrt{1/N}) \] for every \(q=1,\ldots,Q\).
How should the variance be estimated?\ Through a residual technique developed in Deville and Till’e (2005). This technique is comparable to the technique used to calculate the variance of the calibration estimator and has been validated through a set of simulations.
Can balancing and calibration be used simultaneously?\ Both techniques can be used together. There is no contradiction. The best sampling strategy would consist of using them together. In fact, calibration can fix the rounding problem after balancing. Moreover, different variables can be used in calibration from those used in balancing.
What software should be used?\ In SAS-IML, two packages exist (INSEE and University of Neuchatel); in R, the sampling package allows the cube method to be used. This software is available online free of charge.

11.4 Exercises

Suppose a sampling design of size \(n=2\) for a population of size \(N=3\) with an auxiliary information characteristic such that \(x_k=\pi_k\) (k=1,2,3) and, additionally, \(\pi_1+\pi_2+\pi_3=2\).

Write the balancing equations.
Calculate the entries of the matrix \(\mathbf{A}\) (Section 15.2.1).
Define the null space of the matrix \(\mathbf{A}\); that is, \(\ker(\mathbf{A})\).
Obtain the explicit form of \(Q=\boldsymbol{\pi}+\ker(\mathbf{A})\).

Suppose a balanced sampling design with \(N=8\) and \(n=4\). Assume that the first-order inclusion probability vector is \[\boldsymbol{\pi}=\left(\frac{1}{9},\frac{2}{9},\frac{3}{9},\frac{4}{9},\frac{5}{9},\frac{6}{9},\frac{7}{9},\frac{8}{9}\right)'\] and that there are two balancing variables: the first, \(x_{1k}=\pi_k\), and the second, \(x_{2k}=1\), for all \(k\in U\).

Write the balancing equations.
Calculate the entries of the matrix \(\mathbf{A}\).
If the cost function is \[Cost_1(\mathbf{s})=\sum_{p=1}^P\frac{(\hat{t}_{x_p,\pi}-t_x)^2}{t_x^2}\] Obtain the cost generated by the landing phase for the samples:\ \(\mathbf{s}_1=(1,0,0,0,0,1,1,1)'\).\ \(\mathbf{s}_2=(0,0,0,1,1,1,0,1)'\).\ \(\mathbf{s}_3=(0,0,1,1,0,0,1,1)'\).\ \(\mathbf{s}_4=(0,0,1,1,0,1,1,0)'\).\
If the cost function is \[Cost_2(\mathbf{s})=(\mathbf{s}-\boldsymbol{\pi})'\mathbf{A}'(\mathbf{AA}')^{-1}\mathbf{A}(\mathbf{s}-\boldsymbol{\pi})\] Obtain the cost generated by the landing phase for the previous samples.

Prove or refute the following statements:

“Using balanced sampling always improves the efficiency of the sampling strategy.”
“Using calibration always improves the efficiency of the balanced sampling strategy.”
“Using calibration and balanced sampling always improves the efficiency of the sampling strategy.”

Ardilly, P. 1991. “Échantillonnage Représentatif Optimum - Probabilités Inégales.” Annales d’économie Et de Statistique 23: 91–113.

Deville, J. C. 1992. “Constrained Samples , Conditional Inference, Weihting: Three Aspects of the Utilisation of Auxiliary Information.” In Proceedings of the Workshop on the Uses of Auxiliary Information in Survey, edited by Sweeden Örebro.

Deville, J. C., and Y. Tillé. 2004. “Efficient Balanced Sampling: The Cube Method.” Biometrika 91: 893–912.

Deville, J.-C, and Y. Till’e. 2005. “Variance Approximation Under Balanced Sampling.” Journal of Statistical Planning and Inference 128: 411–25.

Neyman, J. 1934. “On the Two Differents Aspects of the Representative Method: The Method of Stratified Sampling and the Method of Purposive Selection.” Journal of the Royal Statistical Society 97: 558–625.

Royal, R. M., and J. Herson. 1973. “Robust Estimation in Finite Population II: Estratification on a Size Variable.” Journal of the American Statistical Association 68: 891–93.

Särndal, C. E., B. Swensson, and J. Wretman. 1992. Model Assisted Survey Sampling. Springer, New York.

Till’e, Y. 2006. Sampling Algorithms. Springer.

Valliant, R., A. H. Dorfman, and R. M. Royall. 2000. Finite Population Sampling and Inference. Wiley.

Yates, F. 1946. “A Review of Recent Statistical Developments in Sampling and Sampling Surveys.” Journal of the Royal Statistical Society A109: 12–43.