# Enhanced Cube Implementation
For Highly Stratified Population

Raphaël Jauslin111Institute of statistics, University of Neuchâtel, Bellevaux
51, 2000 Neuchâtel, Switzerland (E-mail<EMAIL_ADDRESS>, Esther
Eustache11footnotemark: 1 and Yves Tillé11footnotemark: 1

###### Abstract

A balanced sampling design should always be the adopted strategies if
auxiliary information is available. Besides, integrating a stratified
structure of the population in the sampling process can considerably reduce
the variance of the estimators. We propose here a new method to handle the
selection of a balanced sample in a highly stratified population. The method
improves substantially the commonly used sampling design and reduce the time-
consuming problem that could arise if inclusion probabilities within strata do
not sum to an integer.

Key words: balanced sampling, clustered sampling, auxiliary information,
unequal probability sampling

##  1 Introduction

In survey statistics, balanced sampling is a particularly efficient method
when values of auxiliary variables are available for all units in the
population. The idea is to select the sample so that the totals of the
Horvitz-Thompson estimators of some auxiliary variables equal the population
totals. There are different methods for selecting a balanced sample. Deville
and Tillé, (2004) have proposed the cube method which successively transforms
the vector of inclusion probabilities into a sample. The method has been
improved by Chauvet and Tillé, (2006) by reducing the computation time.

In many areas, it is very useful to use stratified sampling designs. As
already indicated by Neyman, (1934), the variance of the Horvitz-Thompson
estimator can be reduced by constructing strata such that the variables are
homogeneous within the strata. Besides, Chauvet, (2009) proposed a specific
algorithm to obtain balanced samples in the strata of a population. However,
this method becomes cumbersome when the number of strata is large.

A highly stratified population is very common in survey sampling. For example,
it may be necessary to select individuals from a population while requiring
that at most only one individual from each household in a population is taken.
Each household is then a stratum. In spatial statistics, one can also
construct small strata of neighboring units to obtain well-spread samples.
Highly stratified sampling is also necessary for some donor imputation
methods: the objective is to select a respondent for each non-respondent to
impute its values. Each non-respondent then defines a stratum from which a
respondent is to be selected (Hasler and Tillé,, 2016).

The balanced and stratified sampling method of Chauvet, (2009) has been
improved by Hasler and Tillé, (2014) to partially resolve the disadvantage of
the time required to process a highly stratified population. When the sum of
the inclusion probabilities in the strata is not an integer, the computation
time can become problematic. This problem arises, for example, when the
objective is to select less than one individual per household. Neither of the
two methods already proposed solves the computational time problem in these
situations.

In this paper, we propose a new method to obtain a stratified balanced sample.
This new method is particularly interesting when the population is highly
stratified and the inclusion probabilities do not sum to an integer within the
strata. We refer readers to Tillé, (2020) and Hankin et al., (2019) to have
more information on the general settings on stratified balanced sampling
design.

The document is organized as follows. The section 2 gives the basic notations
and settings. Section 3 present the problem of selecting a balanced sample. In
the section 4, we review the cube method and how it is used to select a
balanced sample. In section 5, we discuss the issue of the highly stratified
population and review the methods used to select a sample in this case. In the
section 6, we present the new method and the section 7 is devoted to variance
estimation. In the section 8, we give the simulation results of the different
algorithms on an artificial dataset while the section 9 gives a conclusion on
the new method.

##  2 Basic sampling notations

Consider a finite population $U$ of size $N$ whose units can be defined by
labels $k\in\\{1,2,\dots,N\\}$. Let define a variable of interest $y$. Suppose
that we are trying to estimate the following unknown total:

$Y=\sum_{k\in U}y_{k}.$ (1)

A sampling design is defined by the probability $p(s)$ of selecting each
possible subset $s\subset U$ such that $\sum_{s\subset U}p(s)=1.$ Consider a
vector $\textbf{a}=(a_{1},\dots,a_{N})^{\top}$ that maps elements of a subset
$s$ to an $N$ vector of 0s and 1s such that:

$a_{k}=\left\\{\begin{array}[]{lll}1&\text{ if }k\in s,\\\ 0&\text{
otherwise},\end{array}\right.$

for $k\in U$. For each unit of the population, the inclusion probability
$\pi_{k}$, with $0\leq\pi_{k}\leq 1$, is defined as the probability of
selecting $k$ into a sample $s$:

$\pi_{k}=\textrm{P}(k\in s)=\textrm{E}(a_{k})=\sum_{s\subset U|k\in s}p(s).$

Let $\boldsymbol{\pi}=(\pi_{1},\dots,\pi_{N})^{\top}$ be the vector of all the
inclusion probabilities. Let also $\pi_{k\ell}$ be the probability of
selecting units $k$ and $\ell$ together in the sample, with
$\pi_{kk}=\pi_{k}$. Assuming that $\pi_{k}>0$ for all $k\in U$, the total (1)
can be estimated using the classical unbiased Horvitz-Thompson estimator
defined by

$\widehat{Y}=\sum_{k\in U}\frac{y_{k}a_{k}}{\pi_{k}}.$ (2)

##  3 Stratified balanced sampling

Usually, some auxiliary information are available for each unit $k\in U$ in a
vector ${\bf x}_{k}=(x_{k1},\dots,x_{kq})^{\top}\in\mathbb{R}^{q}$, with
$q\in\mathbb{N^{*}}$. A sampling design is said to be balanced on the $q$
auxiliary variables if and only if it satisfies the following balancing
equation:

$\widehat{{\bf X}}=\sum_{k\in s}\frac{{\bf x}_{k}}{\pi_{k}}=\sum_{k\in
U}\frac{{\bf x}_{k}a_{k}}{\pi_{k}}=\sum_{k\in U}{\bf x}_{k}={\bf X}.$

Sometimes, selecting a sample that satisfies exactly the constraints is not
possible due to the rounding problem.

In many applications, inclusion probabilities are such that the selected
sample has a fixed size. In order to obtain a sampling design with fixed
sample size, a linear combination of the auxiliary variables must be
proportional or equal to the vector of inclusion probabilities, i.e. there
exists $\boldsymbol{\psi}\in\mathbb{R}^{q}$ such that
$\boldsymbol{\psi}^{\top}{\bf x}_{k}=\pi_{k}$, for all $k\in U$. Indeed, this
gives

$\sum_{k\in s}\frac{\boldsymbol{\psi}^{\top}{\bf x}_{k}}{\pi_{k}}=\sum_{k\in
s}\frac{\pi_{k}}{\pi_{k}}=n.$

The size of the sample will be fixed only if $n$ is an integer. If it is not
the case, the sample size will be equal to the higher or lower integer to $n$.

More generally, the problem of selecting a balanced sample is written as the
following linear system :

$\left\\{\begin{array}[]{lll}\displaystyle{\bf A}^{\top}{\bf a}={\bf
A}^{\top}\boldsymbol{\pi},\\\ {\bf a}\in\\{0,1\\}^{N},\end{array}\right.$ (3)

where ${\bf A}=\left({\bf x}_{1}/\pi_{1},\dots,{\bf
x}_{N}/\pi_{N}\right)^{\top}$. The aim consists then of obtaining a sample
${\bf a}$ that satisfies (or approximately satisfies) the constraints.

Suppose that the population $U$ is divided into $H$ strata
$U_{1},\dots,U_{H}$, with respective sizes of $N_{1},\dots,N_{H}$. The strata
form a partition and respect the following properties:

$U=\bigcup_{h=1}^{H}U_{h},~{}N_{h}>0,~{}U_{h}\cap U_{\ell}=\emptyset,\text{
for all }h,\ell\in\\{1,\dots,H\\}$

Then, this implies that $\sum_{h=1}^{H}$ $N_{h}=N$. The inclusion
probabilities sum to a value $n_{h}$ in each stratum $h$, i.e.
$n_{h}=\sum_{k\in U_{h}}\pi_{k}$. Let ${\bf h}=(h_{1},$ $\dots,$
$h_{N})^{\top}$ be a categorical vector that specifies the stratum to which
each unit belongs. For example, $h_{k}=\ell$ means that unit $k$ belongs to
strata $U_{\ell}$, with $k\in U$ and $\ell\in\\{1,\dots,H\\}$. Another way for
expressing the stratum of each unit is to use the disjunctive form. Let ${\bf
H}$ be the disjunctive matrix of the corresponding vector ${\bf h}$ of size
$N\times H$, such that:

$\displaystyle{\bf H}=\big{(}{\bf 1}(U_{1}),\dots,{\bf 1}(U_{H})\big{)},$

where ${\bf 1}(U_{h})\in\mathbb{R}^{N}$ is an column vector such that its
$k$th element is equal to 1 if the unit $k$ belongs to the stratum $U_{h}$ and
0 otherwise.

Obtaining a balanced sample in a stratified population is equivalent to adding
stratification constraints to the previous linear system (3). These
constraints are contained in the matrix ${\bf H}$, so the modification of the
linear problem gives:

$\left\\{\begin{array}[]{lll}\displaystyle\left({\bf H}\ {\bf
A}\right)^{\top}{\bf a}=\left({\bf H}\ {\bf
A}\right)^{\top}\boldsymbol{\pi},\\\ {\bf
a}\in\\{0,1\\}^{N}.\end{array}\right.$ (4)

The number of constraints in the linear problem is then $\left(q+H\right)$. In
the next section, a method to select a balanced sample is presented.

##  4 Cube Method

Deville and Tillé, (2004) developed the cube method that selects a balanced
sample respecting the inclusion probabilities. The method can deal with equal
or unequal inclusion probabilities. The algorithm is separated into two
phases.

  * •

The first phase is called the flight phase. It modifies recursively and
randomly the vector of inclusion probabilities $\boldsymbol{\pi}$ into a
sample by respecting exactly the balancing constraints of the problem. The
subspace induced by the linear system (3) could be rewritten using the
following notation:

$\mathcal{A}=\left\\{{\bf a}\in\mathbb{R}^{N}|{\bf A}^{\top}{\bf a}={\bf
A}^{\top}\boldsymbol{\pi}\right\\}=\boldsymbol{\pi}+\text{Null}({\bf
A}^{\top}),$

where $\text{Null}({\bf A}^{\top})=\left\\{{\bf u}\in\mathbb{R}^{N}|{\bf
A}^{\top}{\bf u}=0\right\\}$. The idea is then to use a vector ${\bf u}$ of
the null space of ${\bf A}^{\top}$ in order to update randomly the vector
$\boldsymbol{\pi}$. The whole procedure of the update can be found in Deville
and Tillé, (2004). At each step, at least one component is set to 0 or 1.
Matrix ${\bf A}$ is updated with the new inclusion probabilities. This step is
repeated until the null space of ${\bf A}^{\top}$ is empty. At the end of the
flight phase, the final updated vector of $\boldsymbol{\pi}$ contains at most
$q$ elements that are still not equal to 0 or 1.

  * •

The second phase is called the landing phase. This phase allows to obtain the
sample ${\bf a}$ that respects as much as possible the balancing constraints.
There are two different ways to achieve it, by relaxing the $q$ constraints
one by one, or by linear programming.

In the flight phase, the major computational cost comes from the research of a
vector in the null space of ${\bf A}^{\top}$. Chauvet and Tillé, (2006) have
improved this time-consuming inconvenience using a sub-matrix of ${\bf A}$
rather than the entire matrix. The idea is to consider a submatrix that has
one more row than the number of columns to ensure to have at least one vector
in its null space. This submatrix, denoted by ${\bf B}$, has then a size of
$(q+1)\times q$, with respect to $q<N$ and $\mbox{Rank}({\bf B})\leq q$.

The interest of using this submatrix comes from the following result: a vector
${\bf u}$ of $\text{Null}({\bf B}^{\top})$ completed by $\left(N-(q+1)\right)$
zeros is a vector of $\text{Null}({\bf A}^{\top})$. With this idea, all the
computations can be done using only a submatrix ${\bf B}$. Usually, $N$ is
much greater than $q$, the size of ${\bf B}$ is then much smaller than ${\bf
A}$. This implies obviously an important gain of computational time. The
method proposed in this paper uses the same idea. In the next section, the
particular case of highly stratified sampling is considered.

##  5 Highly stratified population

It is always preferable to consider a stratified population in order to
estimate the total (1). Indeed, the variance of the Horvitz-Thompson estimator
(2) can be considerably reduced compared to the non-stratified estimator (1).
However, when the population is highly stratified (i.e. $H$ is very large),
the selection of a balanced sample with classical methods becomes difficult
due to the too large number of constraints in ${\bf H}$. In order to decrease
the time-consuming problem, different approaches have already been proposed.

Chauvet, (2009) has developed an algorithm to select a balanced sample in a
highly stratified population. Firstly, a flight phase is applied inside each
stratum. This allows modifying the inclusion probabilities such that these are
as balanced as possible in each stratum. Next, a flight phase is applied on
the whole population. Finally, a landing phase is carried out on units that
are not still selected or rejected. This procedure has the advantage to be
simple to implement. Its major deficiency is when the number of strata $H$
becomes too large, the procedure remains very slow and often cannot even be
used.

Hasler and Tillé, (2014) have proposed another method to deal with highly
stratified population. As the previous method, it begins by applying the
flight phase of the cube method to each stratum of the population. Next, it
carries out a flight phase on an union of strata by adding another stratum at
each step. By doing this, strata are managed one after the other and the
inclusion probabilities of certain strata are set to 0 or 1 during this step.
The idea behind this procedure is to reduce the matrix ${\bf H}$ considered
because some strata are removed from the matrix when all its units are
selected or rejected. At the end, a landing phase is applied. However, if
$n_{h}$ is not equal to an integer for a stratum $U_{h}$, this method also
remains very time-consuming. Indeed, some strata are never completely removed
during the procedure and then the submatrix of ${\bf H}$ considered becomes
too large.

The properties of the cube method imply that the inclusion probabilities are
satisfied and that the sample is balanced on the auxiliary variables in these
two methods. However, they still have difficulty to deal with all the
situations of highly stratified sampling. In the next following section, a new
method is presented in order to completely resolve these drawbacks.

##  6 Proposed method

In the fast implementation of the cube method (Chauvet and Tillé,, 2006), the
main modification was to use a matrix smaller than ${\bf A}$ to update
$\boldsymbol{\pi}$. This allows to considerably decreasing the computational
cost. The idea of our method is similar but adapted to a stratified
population: consider a matrix of constraints ${\bf B}$ smaller than $({\bf H}\
{\bf A})$ during the use of the cube method.

The submatrix matrix ${\bf B}$ must be found at each step of the flight phase
of the cube method. As explained in Section 3, the number of balancing
constraints depends on the number of strata $H$ when the population is
stratified. By considering a matrix ${\bf B}$ with fewer rows, or units, the
corresponding vector of strata ${\bf h}$ will be reduced. This subvector of
${\bf h}$ will contain fewer categories and then the corresponding matrix
${\bf H}$ will have fewer columns. The number of constraints will therefore
depend on the rows of ${\bf B}$. This is why obtaining the matrix ${\bf B}$
with exactly one row more than its number of columns is not as easy as with an
unstratified population. Algorithm 1 shows how to find the number of rows to
consider to obtain the smaller matrix ${\bf B}$ such that ${\bf B}$ has
exactly one row more than its number of columns.

Algorithm 1 Find the sub-matrix ${\bf B}$ of $\left({\bf H}\ {\bf A}\right)$

Let $q$ be the number of auxiliary variables of ${\bf A}$. Initialize $q^{1}$
by $q$. For $t=1,2,3,\dots$ repeat the following steps:

  1. 1.

Extract the first $q^{t}$ rows of the vector ${\bf h}$ and denote it ${\bf
h}^{t}$.

  2. 2.

Denote $H^{t}$ the number of different strata in ${\bf h}^{t}$.

  3. 3.

Update $q^{t+1}=q+H^{t}+1$.

while $q^{t+1}>q^{t}$.
Finally, ${\bf B}$ is defined as the $q^{t}$ first rows of the concatenated
matrix $\left({\bf H}^{t}\ {\bf A}^{t}\right)$, where ${\bf A}^{t}$ and ${\bf
H}^{t}$ are the submatrix containing only its $q^{t}$ first rows.

###### Example 6.1

Suppose that $q=2$ and that the categorical vector is equal to ${\bf h}=(1,$
$1,$ $2,$ $2,$ $3,$ $3,$ $3,$ $4,$ $4)^{\top}$. We obtain

$\begin{array}[]{llllll}t=1:&q^{1}=2,&{\bf
h}^{1}=(1,1)^{\top},&H^{1}=1&\to&q^{2}=2+1+1=4,\\\ t=2:&q^{2}=4,&{\bf
h}^{2}=(1,1,2,2)^{\top},&H^{2}=2&\to&q^{3}=2+2+1=5,\\\ t=3:&q^{3}=5,&{\bf
h}^{3}=(1,1,2,2,3)^{\top},&H^{3}=3&\to&q^{4}=2+3+1=6,\\\ t=4:&q^{4}=6,&{\bf
h}^{4}=(1,1,2,2,3,3)^{\top},&H^{4}=3&\to&q^{5}=2+3+1=6,\\\
t=5:&q^{5}=q^{4}\end{array}$

${\bf B}$ contains then $q^{4}=6$ rows and 2 + 3 = 5 columns. So it is a
matrix with only one more rows than the number of columns as desired.

The matrix ${\bf B}$ is found after having computed its number of rows $q^{t}$
using Algorithm 1. The first $q^{t}$ elements of ${\bf h}$ composed the strata
membership vector ${\bf h}^{t}$. The disjunctive matrix ${\bf H}^{t}$ can then
be found using ${\bf h}^{t}$. The matrix ${\bf B}$ is equal to $\left({\bf
H}^{t}\ {\bf A}^{t}\right)$, with ${\bf A}^{t}$ the submatrix of ${\bf A}$
containing only its $q^{t}$ first rows. The same procedure proposed by Chauvet
and Tillé, (2006) can then be applied. If the population is highly stratified
and the number of auxiliary variables is acceptable, our procedure can be very
efficient. Moreover, it handles inclusion probabilities that do not sum to an
integer inside strata. Algorithm 2 presents the whole method.

Algorithm 2

Consider $\boldsymbol{\pi}$ the $N$ vector of inclusion probabilities such
that $0<\pi_{k}<1$, for $k\in\\{1,\dots,N\\}$.

  * I.

Perform a flight phase on each stratum according to the inclusion
probabilities $\boldsymbol{\pi}$ and the balancing constraints in ${\bf
A}^{\top}$. The vector $\boldsymbol{\pi}$ is updated by $\boldsymbol{\pi}^{1}$
such that some of its elements are set to 0 or 1. Compute the set of indices
${\bf i}^{1}\subset\\{1,\dots,N\\}$ containing the unit indices that have an
inclusion probability still not equal to 0 or 1.

  * II.

Initialize $t$ by 1. Repeat step 1. to 6. until it is no more possible to find
the matrix ${\bf B}$ or until the vector ${\bf u}$ is null.

    1. 1.

In ${\bf A}$, ${\bf h}$ and $\boldsymbol{\pi}$, consider only units with
indices in ${\bf i}^{t}$.

    2. 2.

Apply the Algorithm 1 to find the submatrix ${\bf B}$ of $({\bf H}\ {\bf A})$.

    3. 3.

Compute ${\bf u}$, a vector of the null space of ${\bf B}$ completed by 0s to
obtain a vector with the same size as ${\bf i}^{t}$.

    4. 4.

Compute $\lambda_{1}>0$ and $\lambda_{2}>0$, the two greater values such that

$\begin{array}[]{ccccc}0&\leqslant&\pi_{k}^{t}+\lambda_{1}u_{k}\leqslant 1\\\
0&\leqslant&\pi_{k}^{t}-\lambda_{2}u_{k}\leqslant 1\\\ \end{array}\text{, for
all }k\in{\bf i}^{t}.$

    5. 5.

Update $\boldsymbol{\pi}^{t}$ by:

$\boldsymbol{\pi}^{t+1}=\left\\{\begin{array}[]{cccc}\boldsymbol{\pi}^{t}+\lambda_{1}{\bf
u}&\text{ with probability }&\lambda_{2}/(\lambda_{1}+\lambda_{2}),\\\
\boldsymbol{\pi}^{t}-\lambda_{2}{\bf u}&\text{ with probability
}&\lambda_{1}/(\lambda_{1}+\lambda_{2}).\end{array}\right.$

    6. 6.

Update $t$ by $t+1$ and update ${\bf i}^{t}$ the set of indices containing the
unit that have an inclusion probability still not equal to 0 or 1.

  * III.

It could remain some units that are still not rejected or selected. Perform a
landing phase by suppression of variables on the balancing variables $({\bf
H}\ {\bf A})$ on the remaining indices ${\bf i}^{t}$ .

##  7 Variance estimation

The variance can be approximated using the method proposed by Deville and
Tillé, (2005). Let the vector

${\bf z}_{k}=({\bf H}~{}~{}{\bf A})_{k},$

where $({\bf H}~{}~{}{\bf A})_{k}$ denote the $k$th row of the matrix $({\bf
H}~{}~{}{\bf A})$. The variance of the Horvitz-Thompson estimator of the total
$\widehat{Y}$ can be approximated by

$\text{var}_{app}(\widehat{Y})=\sum_{k\in
U}c_{k}\left(\dfrac{y_{k}}{\pi_{k}}-\boldsymbol{\alpha}^{\top}{\bf
z}_{k}\right)^{2},$ (5)

where

$c_{k}=\pi_{k}(1-\pi_{k})\frac{N}{N-(H+q)}\text{ and
}\boldsymbol{\alpha}=\left(\sum_{\ell\in U}c_{\ell}{\bf z}_{\ell}{\bf
z}_{\ell}^{\top}\right)^{-1}\sum_{\ell\in U}c_{\ell}{\bf
z}_{\ell}\dfrac{y_{\ell}^{\top}}{\pi_{\ell}}.$ (6)

There exists many different ways to express the quantity $c_{k}$ and then this
leads to various approximations of the variance. Value $c_{k}$ can in
particular be approximated by

$\widetilde{c_{k}}=(1-\pi_{k})\frac{n}{n-(H+q)}.$

Equation (5) can be estimated on a sample $s$ using $\widetilde{c_{k}}$
instead of $c_{k}$ and by replacing the sums on $U$ by sums on $s$ in
Expressions (5) and (6).

##  8 Simulations

In this section, the performance of the method is evaluated on real data
produced by the Swiss Federal Statistical Office, (2020). The dataset contains
information on Swiss establishments. We restrict the study to the Switzerland
region called Espace Mittelland (a region of the second degree of the
Nomenclature of Territorial Units for Statistics (NUTS) of Switzerland). This
region contains 5 cantons (a region of the third degree of the NUTS) and 675
municipalities. For confidentiality reasons, the units considered are the
hectares of land in which at least one establishment is located. In order to
be able to estimate the variance, only 3 hectares of land per municipalities
are included in the study. This implies that the dataset contains information
from 2025 hectares including at least one establishment.

We stratify the units in two different ways: by cantons and by municipalities.
The number of strata is then respectively equal to $H_{c}=5$ and $H_{m}=675$.
Figure 1 shows the dataset with the two proposed stratification. The idea
behind this procedure is to compare the execution time for a stratified
population with a low number of strata versus a high one.

$46.4^{\circ}$N$46.6^{\circ}$N$46.8^{\circ}$N$47.0^{\circ}$N$47.2^{\circ}$N$47.4^{\circ}$N$6.5^{\circ}$E$7.0^{\circ}$E$7.5^{\circ}$E$8.0^{\circ}$E$8.5^{\circ}$ELongitudeLatitude$N_{h}$931593274081038Cantons$46.4^{\circ}$N$46.6^{\circ}$N$46.8^{\circ}$N$47.0^{\circ}$N$47.2^{\circ}$N$47.4^{\circ}$N$6.5^{\circ}$E$7.0^{\circ}$E$7.5^{\circ}$E$8.0^{\circ}$E$8.5^{\circ}$ELongitudeLatitude3$N_{h}$Municipalities
Figure 1: Data extract from the Swiss establishement data base of the Swiss
Federal Statistical Office Swiss Federal Statistical Office, (2020). The data
are restricted to the NUTS region 2. The left plot is showing the separation
by Canton $H_{c}=5$, the right one the separation by Municipality $H_{m}=675$.
The grey gradient scale give the number of units considered in each Canton.
The data are selected such that each municipality contain 3 units.

To compare the method, we will use balancing variables ${\bf x}_{j}$
containing the number of women employed in a sector $j$, with $j=1,2,3$. Each
sector represents a type of activity: sector $j=1$ involves the natural
environment and agriculture, sector $j=2$ is manufacture and sector $j=3$ is
related to services. Table 1 shows the mean time execution of three methods
for highly stratified sampling: the methods of Hasler and Tillé, (2014) and
Chauvet, (2009), detailed in Section 5, and the new one presented in this
article.

Inside each stratum, the inclusion probabilities are equal. For each
stratification, we carry out two different sampling: one with inclusion
probabilities that sum to an integer number within each stratum and one with a
non-integer sum. Then, we consider $n_{h}=80$ and $n_{h}=80.4$,
$h\in\\{1,\dots,H_{c}\\}$, for the first stratification. For the second one,
we take $n_{h}=2$ and $n_{h}=1.4$, $h\in\\{1,\dots,H_{m}\\}$. We choose to
deal also with non-integers $n_{h}$ with the aim to compare the impact of this
situation on the mean sampling time.

Chauvet’s method cannot be compared because its execution time is too long and
should be avoided for highly stratified population. However, it remains very
efficient if the number of strata is acceptable. If inclusion probabilities in
stratum sum to an integer, the Hasler’s method performs very well. However,
the execution time increases strongly when $n_{h}$ is not integer. The
proposed method has a very well behaved and the time is considerably reduced
for a highly stratified population.

In order to compare the variance of the method with the others, we estimate
the variance using some variables of interest $y_{j}$ that contains the total
number of employee of the sector $j$, $j=1,2,3$. In Table 1, we compare the
approximated variance (5), the estimated variance and the simulated variance
computing using the equation:

$v_{sim}=\frac{1}{m}\sum_{s}\\{\widehat{Y}(s)-Y\\}^{2},$ (7)

where $m$ is the number of simulations.

For each method, we vary the number of selected units within each stratum by
taking $n_{h}$ equal to 2 for the stratification with $H_{c}$ strata and 80
for the stratification with $H_{m}$ strata. This implies sample of size
respectively equal to 400 and 1350. The variance estimator seems to be
unbiased for the approximated variance. However, we see that the approximated
variance and estimator are slightly biased to the $v_{i}$. This coming from
the landing phase of each method. We can conclude that the proposed method is
comparable in terms of variance to other methods.

Table 1: Results of 1000 simulations on the Swiss establishments dataset. The
population of size is equal to 2025. We compute the mean time execution in
secondes of each sampling procedure. We vary the number of strata $H$ and the
number of unit selected within each strata $n_{h}$.

| Algorithm
---|---
| Proposed method | Hasler’s method | Chauvet’s method
Cantons ($H=5$)
$n_{h}=80$ | 0.24 | 0.25 | 0.24
$n_{h}=80.4$ | 0.24 | 0.24 | 0.24
Municipalities ($H=675$)
$n_{h}=2$ | 0.42 | 0.4 | 418.07
$n_{h}=1.4$ | 0.53 | 400.55 | 701.77

Table 2: Results of 1000 simulations on a population of size is equal to 2025.
The number of strata is equal to 5 for Cantons and 675 for the Municipalities.
For each variable of interest $y_{j}$, $j=1,2,3$ and for each sampling
methods, we compute the variance approximated by the simulations (7) ,
approximated variance (5) as well as the variance estimator (6)

| Algorithm
---|---
| Proposed method | Hasle’s method | Chauvet’s method
| $v_{sim}$ | $\widehat{\text{var}}(\widehat{Y})$ | $\text{var}_{app}(\widehat{Y})$ | $v_{sim}$ | $\widehat{\text{var}}(\widehat{Y})$ | $\text{var}_{app}(\widehat{Y})$ | $v_{sim}$ | $\widehat{\text{var}}(\widehat{Y})$ | $\text{var}_{app}(\widehat{Y})$
Cantons $(H=5)$
$y_{1}$ | 25788.227 | 24365.755 | 24735.596 | 24286.279 | 24147.542 | 24735.596 | 22913.046 | 24368.969 | 24735.596
$y_{2}$ | 136146.545 | 123257.656 | 131861.315 | 130340.065 | 121443.551 | 131861.315 | 134082.01 | 124254.374 | 131861.315
$y_{3}$ | 145849.824 | 129921.58 | 135199.514 | 148338.808 | 128108.505 | 135199.514 | 153747.141 | 129142.588 | 135199.514
Municipalities $(H=675)$
$y_{1}$ | 1656.128 | 1726.496 | 1736.933 | 1815.448 | 1727.489 | 1736.933 | 1771.027 | 1731.304 | 1736.933
$y_{2}$ | 9177.576 | 9394.483 | 9435.185 | 9109.812 | 9427.651 | 9435.185 | 9887.571 | 9408.017 | 9435.185
$y_{3}$ | 9033.066 | 8990.18 | 8982.537 | 9589.876 | 8950.333 | 8982.537 | 10533.724 | 8937.991 | 8982.537

##  9 Conclusion

The stratified sampling procedure is a well-known and appropriate procedure to
reduce the variance of the Horvitz-Thompson estimator. In this paper, we
propose a new method and implementation that provides an excellent executive
time and flexibility that the existing methods did not allow. In many surveys
where the population is stratified, the sum of inclusion probabilities within
each stratum is not an integer. Other methods are not directly applicable in
this case. We have shown by mean of simulations that the variance of the
estimator is not impacted by our method. All of these results indicate that
our proposed algorithm is very efficient to select a sample in a stratified
and highly stratified population.

## References

  * Chauvet, (2009) Chauvet, G. (2009).  Stratified balanced sampling.  Survey Methodology, 35:115–119.
  * Chauvet and Tillé, (2006) Chauvet, G. and Tillé, Y. (2006).  A fast algorithm of balanced sampling.  Journal of Computational Statistics, 21:9–31.
  * Deville and Tillé, (2004) Deville, J.-C. and Tillé, Y. (2004).  Efficient balanced sampling: The cube method.  Biometrika, 91:893–912.
  * Deville and Tillé, (2005) Deville, J.-C. and Tillé, Y. (2005).  Variance approximation under balanced sampling.  Journal of Statistical Planning and Inference, 128:569–591.
  * Hankin et al., (2019) Hankin, D., Mohr, M., and Newman, K. (2019).  Sampling Theory: For the Ecological and Natural Resource Sciences.  Oxford University Press, New York.
  * Hasler and Tillé, (2014) Hasler, C. and Tillé, Y. (2014).  Fast balanced sampling for highly stratified population.  Computational Statistics and Data Analysis, 74:81–94.
  * Hasler and Tillé, (2016) Hasler, C. and Tillé, Y. (2016).  Balanced $k$-nearest neighbor imputation.  Statistics, 105:11–23.
  * Neyman, (1934) Neyman, J. (1934).  On the two different aspects of the representative method: The method of stratified sampling and the method of purposive selection.  Journal of the Royal Statistical Society, 97:558–606.
  * Swiss Federal Statistical Office, (2020) Swiss Federal Statistical Office (2020).  Statistique structurelle des entreprises (STATENT) Description des données GEOSTAT.  Neuchâtel, Switzerland.
  * Tillé, (2020) Tillé, Y. (2020).  Sampling and estimation from finite populations.  Wiley, New York.