# Sequential Spatially Balanced Sampling

Raphaël Jauslin, Bardia Panahbehagh and Yves Tillé111Raphaël Jauslin and Yves
Tillé, Institute of statistics, University of Neuchatel, Av. de Bellevaux 51,
2000 Neuchatel, Switzerland, Email<EMAIL_ADDRESS>Bardia
Panahbehagh, Department of Mathematics, Kharazmi University, Tehran, Iran.

###### Abstract

Sequential sampling occurs when the entire population is not known in advance
and data are obtained one at a time or in groups of units. This manuscript
proposes a new algorithm to sequentially select a balanced sample. The
algorithm respects equal and unequal inclusion probabilities. The method can
also be used to select a spatially balanced sample if the population of
interest contains spatial coordinates. A simulation study is proposed on a
dataset of Swiss municipalities. The results show that the proposed method
outperforms other methods.

Keywords inclusion probability, spread sampling, stream sampling, survey
methods

##  1 Introduction

A sample is said to be balanced if the Horvitz-Thompson estimators of the
totals calculated from the sample are equal or nearly equal to the population
totals. The use of balancing samples was first proposed by Gini and Galvani,
(1929). They selected a sample of 29 districts out of 214 to replicate certain
population means. The selection of the sample was not random and the method
was strongly criticized by Jerzy Neyman. It is now well known that a sample
can be selected randomly and balanced simultaneously. The cube method, which
randomly selects a balanced sample, has been proposed by Deville and Tillé,
(2004). Recently, Leuenberger et al., (2021) proposed an improvement. They
showed that changing the order of the units before running the algorithm can
significantly increase the quality of the sample balance.

In environmental studies particularly, data contains spatial coordinates. When
the data are spatially autocorrelated, it is often more accurate to spread the
sample in space. Moreover, well-spread samples in space are balanced on
auxiliary variables (see Grafström and Lundström,, 2013) even if the targets
parameters are nonlinear in the auxiliary variables. Many sampling methods are
currently used to select a well-spread sample. A well-known one is the
generalized random tessellation sampling proposed by Stevens Jr. and Olsen,
(2004). It maps a multi-dimensional space into a real line and uses one
dimension systematic sampling to select a sample. Grafström et al., (2012)
proposed the local pivotal method that uses repulsion between two nearby units
to ensure that the nearby units are not both selected. More recently, Jauslin
and Tillé, (2020) have proposed adaptation of the cube method using a
contiguity matrix to select an almost perfectly well-spread sample.

Therefore, balanced samples and well-spread samples are strongly related.
According to auxiliary variables, balancing will result of a well-spread
sample. But some methods can even do better. A combination of the cube method
and the local pivotal method (Grafström and Tillé,, 2013) allows the selection
of a sample simultaneously well-spread and balanced. They showed, moreover,
that this method diminishes drastically the variance of the Horvitz-Thompson
estimator. Vallée et al., (2015) present besides an application of the local
cube method for forest inventories. Throughout this manuscript, we write
“spatially balanced sample” or “well-spread sample” to refer to the same
interpretation.

All these methods apply to finite populations. This means that we must have
access to the entire population before using the sampling algorithm. We speak
of streaming data when the data arrive one by one or in groups of units.
Examples are network data from the Internet, environmental or biological
studies over a long period of time. In this case, no algorithm can be used to
select a balanced sample during the experimentation process.

In this manuscript, we propose a new method to select a balanced sample and a
spatially balanced if spatial coordinates are available. Moreover, the method
is sequential, i.e., the algorithm does not need to get access to the entire
population to run. In section 2, we introduce the basic concept of survey
sampling theory. In section 3, we expose the balancing equations formally and
gives insight interest of selecting a sample at the same time balanced and
well-spread. In section 4, we present the outline of the method while in
section 5, we introduce spread measure and how the method select a spatially
balanced sample. Section 6 contains simulation results on a Swiss
municipalities dataset.

##  2 Basic Concepts

Consider a finite population $U=\\{1,2,\dots,N\\}$. Let
$\mathcal{S}=\\{s|s\subset U\\}$ denote the power set of $U$. A sampling
design is defined by a probability distribution $p(.)$ on $\mathcal{S}$ such
that

$p(s)\geq 0\text{ for all }s\in\mathcal{S}\text{ and
}\sum_{s\in\mathcal{S}}p(s)=1.$

A random sample $S$ is a random vector that maps elements of $\mathcal{S}$ to
an $N$ vector of 0 or 1 such that $\textrm{P}(S=s)=p(s)$. Let define
$a_{k}(S)$, for $k=1,\dots,N$:

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

For ease of notation, $S$ is generally implicit and we write $a_{k}(S)=a_{k}$.
Then a sample can be denoted by means of vector
$\textbf{a}^{\top}=(a_{1},a_{2},\dots,a_{N}).$ For each unit of the
population, the inclusion probability $0\leq\pi_{k}\leq 1$ is defined as the
probability that unit $k$ is selected into sample $S$:

$\pi_{k}=\textrm{P}(k\in S)={\rm E}_{a}(a_{k})=\sum_{s\in S|k\in s}p(s),\text{
for all }k\in U,$

where ${\rm E}_{a}(.)$ denotes the expectation under the sampling design.

Let $\boldsymbol{\pi}=(\pi_{1},\dots,\pi_{N})^{\top}$ be the vector of
inclusion probabilities. Then, ${\rm E}_{a}({\textbf{a}})=\boldsymbol{\pi}.$
In many applications, inclusion probabilities are such that samples have a
fixed size $n$. Let the set of all samples that have fixed size equal to $n$
be defined by

$\mathcal{S}_{n}=\left\\{\textbf{a}\in\\{0,1\\}^{N}~{}~{}\bigg{|}~{}~{}\sum_{k=1}^{N}a_{k}=n\right\\}.$

The sample is generally selected with the aim of estimating some population
parameters. Let $y_{k}$ denote a real number associated with unit $k\in U$,
usually called the variable of interest. For example, the total

$Y=\sum_{k\in U}y_{k}$

can be estimated by using the classical Horvitz-Thompson estimator of the
total defined by

$\widehat{Y}_{HT}=\sum_{k\in S}\frac{y_{k}}{\pi_{k}}=\sum_{k\in
U}\frac{y_{k}a_{k}}{\pi_{k}}.$ (1)

##  3 Balanced Sampling

Let $\textbf{x}_{k}=(x_{k1},x_{k2},\dots,x_{kp})^{\top}\in\mathbb{R}^{p}$ be
the auxiliary information available for the unit $k$. A sample $S$ is said to
be balanced of variables $x$, if it satisfies the balancing equations given by

$\sum_{k\in S}\frac{\textbf{x}_{k}}{\pi_{k}}=\sum_{k\in U}\textbf{x}_{k}.$

Let ${\bf
A}=(\textbf{x}_{1}/\pi_{1},\textbf{x}_{2}/\pi_{2},\dots,\textbf{x}_{N}/\pi_{N})^{\top}$
denote the auxiliary variables expanded by the inclusion probabilities. The
selection of a balanced sample can be written as the generation of a random
vector a that satisfies the following linear constrants:

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

Grafström and Tillé, (2013) showed that if the population of interest has
spatial coordinates, the units are generally spatially correlated. It is then
interesting to select a well-spread sample. Let suppose a general linear model
with spatial correlation:

$y_{k}=\textbf{x}_{k}^{\top}\boldsymbol{\beta}+\varepsilon_{k},\text{ for all
}k\in U,$

where $\varepsilon_{k}$ is a random variable that satisfies ${\rm
E}_{M}(\varepsilon_{k})=0$ and
$\text{var}_{M}(\varepsilon_{k})=\sigma^{2}(\textbf{x}_{k})=\sigma_{k}^{2}$,
with $\sigma^{2}(.)$ a Lipschitz continuous function. Spatial correlation is
modelled by the function

$\text{cov}_{M}(\varepsilon_{k},\varepsilon_{\ell})=\sigma_{k}\sigma_{\ell}\rho_{k\ell},\text{
with }k\neq\ell\in U$

where $\rho_{k\ell}$ is a function that decreases when the distance between
$k$ and $\ell$ increases. Operators ${\rm E}_{M}(.),\text{var}_{M}(.)$ and
$\text{cov}_{M}(.)$ are respectively the expectation, the variance and the
covariance under the model. Under this model it is possible to show that the
anticipated variance of the Horvitz-Thompson estimator is

${\rm E}_{p}{\rm E}_{M}(\widehat{Y}_{HT}-Y)^{2}={\rm
E}_{p}\left\\{\left(\sum_{k\in S}\frac{\textbf{x}_{k}}{\pi_{k}}-\sum_{k\in
U}\textbf{x}_{k}\right)^{\top}\boldsymbol{\beta}\right\\}^{2}+\sum_{k\in
U}\sum_{\ell\in
U}\sigma_{k}\sigma_{\ell}\rho_{k\ell}\frac{\pi_{k\ell}-\pi_{k}\pi_{\ell}}{\pi_{k}\pi_{\ell}},$
(2)

where $\pi_{k\ell}={\rm E}(a_{k}a_{\ell})$ are the joint inclusion
probabilities. From the first term of Equation (2), the reduction of the
variance is done by selecting a balanced sample. Whereas, the second term is
reduced if the inclusion probabilities $\pi_{k\ell}$ is small while the
$\rho_{k\ell}$ is large. This means that a sample must be selected in a spread
and balanced manner to minimize the anticipated variance.

Deville and Tillé, (2004) have proposed the cube method to select a balanced
sample. The method first consists in performing a random walk inside the
intersection of the hypercube $[0,1]^{N}$ and the subspace spanned by the
balancing constraints to approach a balanced sample. Secondly a sample is
selected such that the balancing equations are approximately satisfied. The
next section discusses the proposed method. It selects a balanced sample
sequentially, which means that the method could be applied on a stream, the
full population is not needed.

##  4 Outline of the Proposed Method

In this section, we present the main idea of the method. Suppose a population
$U$ where units arrive sequentially, by groups of units or one by one. We want
to respect balancing equations and take a decision for the current unit
observed in the stream. For simplicity, let suppose that the current unit is
the first one. The objective of the algorithm is to update only $J<<N$ units
to compensate the decision taken on the current unit. More formally, we try to
find $v_{k}$, $k=2,\dots,J$, characterized by the following update:

$\left\\{\begin{array}[]{lll}\pi_{1}^{1}=0,&\displaystyle\pi_{k}^{1}=\pi_{k}+v_{k},k=2,\dots,J,\\\
\pi_{1}^{2}=1,&\displaystyle\pi_{k}^{2}=\pi_{k}-\frac{1-\pi_{1}}{\pi_{1}}v_{k},\end{array}\right.$

where $\pi_{k}^{1}$ are the updated inclusion probabilities if the decision is
to omit the current unit $k$, in this case $k=1$, same for $\pi_{k}^{2}$ if we
select the current unit $k$. The balancing equations need to be satisfied, we
obtain the following equality:

$\sum_{k=1}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}\pi_{k}^{1}=\sum_{k=1}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}\pi_{k}^{2}=\sum_{k=1}^{J}\textbf{x}_{k},$

that implies

$\sum_{k=1}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}\pi_{k}^{1}=0+\sum_{k=2}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}(\pi_{k}+v_{k})=\frac{\textbf{x}_{1}}{\pi_{1}}+\sum_{k=2}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}\left(\pi_{k}-\frac{1-\pi_{1}}{\pi_{1}}v_{k}\right)=\sum_{k=1}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}\pi_{k},$

and

$\sum_{k=2}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}v_{k}=\frac{\textbf{x}_{1}}{\pi_{1}}-\frac{1-\pi_{1}}{\pi_{1}}\sum_{k=1}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}v_{k}=\textbf{x}_{1},$

and thus

$\sum_{k=2}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}v_{k}=\textbf{x}_{1}.$

To ensure that the updated inclusion probabilities remain between 0 and 1, we
must also have

$0\leq\pi_{k}^{1},\pi_{k}^{2}\leq 1,\text{ for }k=2,\dots,J,$

which implies the two constraints

$\max\left\\{-\pi_{k},(\pi_{k}-1)\frac{\pi_{1}}{1-\pi_{1}}\right\\}\leq
v_{k}\leq\min\left(1-\pi_{k},\pi_{k}\frac{\pi_{1}}{1-\pi_{1}}\right),k=2,\dots,J.$
(3)

In order to find $v_{k}$, we propose to solve the following program:

$\left\\{\begin{array}[]{ll}\displaystyle\text{maximize}&\displaystyle\sum_{k=2}^{J}v_{k}c_{k}\\\
\text{subject
to}&\displaystyle\sum_{k=2}^{J}\frac{\textbf{x}_{k}}{\pi_{k}}v_{k}=\textbf{x}_{1},\\\
&v_{k}\geq-\min\left\\{\pi_{k},(1-\pi_{k})\frac{\pi_{1}}{1-\pi_{1}}\right\\},\\\
&v_{k}\leq\min\left(1-\pi_{k},\pi_{k}\frac{\pi_{1}}{1-\pi_{1}}\right),\end{array}\right.$
(4)

where $c_{k}$ is a cost function that is supposed to penalize more units that
are far to current one. A naive function is $c_{k}=(J-k)$ which decreases as
$k$ increases. If there is no solution the value of $J$ is incremented by one
and we recompute until a solution is found. The initial value for the variable
$J$ must be larger than the dimension of $\textbf{x}_{k}$. If at least one
auxiliary variable is equal to inclusion probabilities, i.e.
$\textbf{x}_{k}=\pi_{k}$, then the method respects fixed sample size.

Suppose now that some steps have been completed and that the current step is
the $i$th one. Let denote $P_{i}=\\{\pi_{i},\dots,\pi_{i+M}\\}$ the pool of
available units at step $i$ where $M\in\mathbb{N}$. This can be either the
entire population or a small part of it if we are in a stream. The algorithm
will take a unit in $P_{i}$ and take a decision on a certain amount of units
$J$, with $J\leq M$. The constraint on the $v_{k}$ depends on the $\pi_{i}$,
i.e. the inclusion probability of the current unit might have a serious impact
on the potential value that could take $v_{k}$. Figure 1 shows possible
$v_{k}$ against $\pi_{k}$, for different value of $\pi_{i}$. The value
$\pi_{k}$ vary from 0 to 1 and the lower and upper bound of the value $v_{k}$
is plotted. It shows clearly that when $\pi_{i}$ is very small, it ends with a
tight interval for $v_{k}$. On the other hand, as the inclusion probabilities
$\pi_{i}$ increases, we obtain a larger range for $v_{k}$. If inclusion
probability of the current unit $i$ is very small, we might have a sharp
increase of the value $J$. This comes from the fact that the balancing
equations are not symmetric for unequal probability sampling. If a balanced
sample is selected with unequal inclusion probabilities, the complementary
sample is not a balanced sample.

To optimize the problem of sharp increase of the value $J$, we propose to
always select as considered unit the one that has the largest inclusion
probability in the pool $P_{i}$ of available units. By taking the one that has
the largest inclusion probability and reordering with respect with the latter,
we optimize the size of $J$ (See Algorithm 1 to have all details of the
implementation). The method is comparable to the cube method. Depending on the
auxiliary variables, it might possible that the method does not end up
directly with a sample. Some units might remain with inclusion probabilities
not equal to an integer. In this case, we propose to land the method by
launching the same process proposed in the cube method, either by suppression
of variables or by linear programming.

Figure 1: For six different inclusion probabilities $\pi_{i}$, we let vary
$\pi_{k}$ from 0 to 1. On the $y$-axis, the bounds of $v_{k}$ in Equations (3)
are calculated. The coloured area represents the eligible value of $v_{k}$,
the shaded line is the upper bound while the bottom line is the lower bound.

##  5 An Enhancement for Spatial Sampling

In addition to the auxiliary variables $\textbf{x}_{k}\in\mathbb{R}^{p}$,
$k\in U$, let suppose that population contains spatial coordinates
$\textbf{z}_{k}\in\mathbb{R}^{q}$. Population with spatial coordinates are
often spatially correlated. In that case, selecting a well-spread sample is a
common practice (Grafström and Lundström,, 2013). A question arises naturally
in spatial sampling: is there measures to see whether a sample is well-spread?
Stevens Jr. and Olsen, (2003) proposed an index based on the Voronoï polygon.
Let $b_{i}$, $i\in S$, be the sum of inclusion probabilities within the $i$th
Voronoï polygons. They showed that the expected value of $b_{i}$ is equal to 1
and proposed to measure how vary these sums using the following quantity:

$B=\frac{1}{n}\sum_{i\in S}(b_{i}-1)^{2}.$

Moran, (1950) proposes to measure spatial correlations in particular case
where the spatial coordinates are placed on a grid. Tillé et al., (2018);
Jauslin and Tillé, (2020) propose a normalized version that uses the notion of
a contiguity matrix ${\bf W}$. Readers can find more information on how the
spatial weights matrix is computed in Jauslin and Tillé, (2020). Moran $I$
index is given by

$I=\frac{(\textbf{a}-\bar{\textbf{a}}_{w})^{\top}{\bf
W}(\textbf{a}-\bar{\textbf{a}}_{w})}{\sqrt{(\textbf{a}-\bar{\textbf{a}}_{w})^{\top}{\bf
D}(\textbf{a}-\bar{\textbf{a}}_{w})(\textbf{a}-\bar{\textbf{a}}_{w})^{\top}{\bf
G}(\textbf{a}-\bar{\textbf{a}}_{w})}},$ (5)

where a is the sample and

$\bar{\textbf{a}}_{w}=\frac{\textbf{a}^{\top}{\bf W}{\bf 1}}{{\bf
1}^{\top}{\bf W}{\bf 1}},$

${\bf D}$ is the $N\times N$ diagonal matrix containing for each $k$,
$\sum_{\ell\in U}w_{k\ell}$ on its diagonal,

${\bf G}={\bf C}^{\top}{\bf D}{\bf C},~{}~{}~{}~{}{\bf C}={\bf D}^{-1}{\bf
W}-\frac{{\bf 1}{\bf 1}^{\top}{\bf W}}{{\bf 1}^{\top}{\bf W}{\bf 1}},$

and ${\bf 1}$ is a column vector of $N$ ones. These two measures are the ones
that we will use in the Section 6.

The proposed method could be used in an interesting way to select a well-
spread sample. At each step before performing the linear program (4), we
select the unit that have the maximum inclusion probability in the pool
$\pi_{(i)}=\max(\pi_{i},\dots,\pi_{i+M})$. We define the reordered pool
$P_{(i)}=\\{\pi_{(i)},\dots,\pi_{(i+M)}\\}$, such that $\pi_{(i+j)}$ is the
$j$th closest unit to $\pi_{(i)}$. Pool is so reordered with respect to the
distance to the unit that have maximum inclusion probability. Doing that will
then modify primarily inclusion probabilities of units closer to the current
one. If a window has a sum of inclusion probabilities closer to 1 and if the
decision of the $i$th unit is to select it (i.e. inclusion probability is
transformed to 1), it will ensure that nearest neighbours of the $i$th unit
will not be selected. This leads to a selection process that spreads the
sample very well in the considered space. If the population units come from a
stream, it is possible that some units have not yet appeared and that their
spatial coordinates are in the neighbourhood of the $i$th unit. This side
effect cannot be controlled, because, if we do not control which unit appears
in the stream, it is not possible to know whether a future unit might appear
within the considered window.

Algorithm 1 Algorithm of stream balanced sampling

Let $\boldsymbol{\pi}$ be the inclusion probability vector,
$\textbf{z}_{k}\in\mathbb{R}^{q}$ the spatial coordinates and
$\textbf{x}_{k}\in\mathbb{R}^{p}$ the value of the auxiliary variables of the
$k$th unit. For $i=1,2,\dots$, suppose a pool of units
$P_{i}=\\{\pi_{i},\dots,\pi_{i+M}\\}$.

  1. 1.

Find $\pi_{(i)}\in P_{i}$ the maximum inclusion probability value in the pool.
Define $P_{(i)}=\\{\pi_{(i)},\dots,\pi_{(i+M)}\\}$ the pool of units reordered
with respect to the distance of $\pi_{(i)}$ calculated using the spatial
coordinates $\textbf{z}_{k}$, i.e. $\pi_{(i+j)}$ is the $j$th closest unit to
the unit $\pi_{(i)}$.

  2. 2.

Find $J$ and $v_{k},~{}k=i+1,\dots,J$ such that the following linear program
has a solution:

$\left\\{\begin{array}[]{ll}\displaystyle\text{maximize}&\displaystyle\sum_{k=i+1}^{J}v_{k}(J-(k-i))\\\
\text{subject
to}&\displaystyle\sum_{k=i+1}^{J}\frac{\textbf{x}_{(k)}}{\pi_{(k)}}v_{k}=\textbf{x}_{(i)}\\\
&v_{k}\geq-\min\left\\{\pi_{(k)},(1-\pi_{(k)})\frac{\pi_{(i)}}{1-\pi_{(i)}}\right\\},\\\
&v_{k}\leq\min\left(1-\pi_{(k)},\pi_{(k)}\frac{\pi_{(i)}}{1-\pi_{(i)}}\right)\end{array}\right.$

At each step $J$ increases and units are reordered with respect to the
distance of the current unit.

  3. 3.

Inclusion probabilities are modified on $P_{(i)}$. For $k=i,\dots,J$, let
$q=\pi_{(i)}$,

$\left\\{\begin{array}[]{llll}\pi_{(i)}=0&\text{and}&\pi_{(k)}=\pi_{(k)}+v_{k}&\text{with
probability }1-q,\\\
\pi_{(i)}=1&\text{and}&\displaystyle\pi_{(k)}=\pi_{(k)}-\frac{1-\pi_{(i)}}{\pi_{(i)}}v_{k}&\text{with
probability }q,\end{array}\right.$

The decision is then taken for unit $\pi_{(i)}$.

  4. 4.

The current unit $\pi_{(i)}$ and units that have inclusion probabilities
transformed to an integer are removed from the pool and new units are
included.

  5. 5.

Repeat steps 1-4 until it is no longer possible to find a solution to the
linear program or if all inclusion probabilities are modified to an integer
value.

  6. 6.

In general, it remains inclusion probabilities not equal to 0 or 1. A landing
phase is launched by using either the doubly balanced sampling algorithm or
landing phase of the cube method if there are no spatial coordinates.

##  6 Application to Swiss Municipalities

In this section, an analysis is performed on the Swiss municipalities summary
dataset (Swiss Federal Statistical Office,, 2021). The dataset contains
different variables giving information on the 2172 municipalities in
Switzerland. Table 1 shows which variables are selected to perform the
analysis. Two framework are compared. Firstly, a framework that uses equal
inclusion probabilities. Secondly, a framework with unequal inclusion
probabilities set up proportional to the population total in each
municipality. Figure 2 shows the auxiliary variables and display the
correlation. In both design, we analyse two variables of interest, the log
number of establishments in the tertiary sector and the foreigners percentage
within each municipality. These two variables are interesting because on both
there exists a strong spatial correlation. Moreover, the tertiary variable is
also strongly correlated to the auxiliary variables selected for the unequal
design. Figure 3 shows the two variable of interest.

To evaluate the performance of the propose method, we run three other
algorithms for each design. The Horvitz-Thompson estimator (1) is calculated
for each sampling design. Local cube (Grafström and Tillé,, 2013) may be the
most important one. It also selects a spatially balanced and a balanced sample
at the same time and the performance in terms of variance should be at least
comparable with this one. Next, we compute the well-spread sampling method,
the local pivotal method (Grafström et al.,, 2012). The idea is to compare the
performance with a purely spatially balanced method. Finally, we compute the
conditional Poisson sampling of fixed size, also called the max-entropy design
(Hájek,, 1981; Deville and Tillé,, 2000; Chen et al.,, 1994). This sampling
design is the most random selective method. It means that this method can be
used as a reference sampling design to see whether the estimation is improved
by balancing and/or spatially balancing the sample.

Table 2 shows results computed on 10’000 simulations. The variance simulated
is computed using the following formula:

$v_{sim}=\frac{1}{10^{\prime}000}\sum_{m=1}^{10^{\prime}000}\left\\{\widehat{Y}_{HT}(s_{m})-Y\right\\}^{2},$

We observe that the proposed method is better to select a spatially balanced
sample for unequal and equal inclusion probabilities in all cases. Moreover,
the variance is always either reduces or comparable. As expected, the local
pivotal method selects a very well-spread sample with unequal inclusion
probabilities but the variance is higher due to the lack of balancing
equations. The proposed method is implemented in the R package Jauslin et al.,
(2021), while the local cube and the local pivotal could be find in the R
package Grafström and Lisic, (2019).

Table 1: Selected variables on the Swiss municipalities dataset. Variables of
interest
---
Foreigner | Number of foreigners.
Tertiary | Logarithm of the number of establishments in the tertiary sector.
Auxiliary variables equal design
Population | Logarithm of the number of people.
Age 0-19 | Number of people between 0 and 19 years old.
Age 20-64 | Number of people between 20 and 64 years old.
Age > 64 | Number of people older than 64 years old.
Deathrate | Death rate.
Establishment | Number of establishments.
Table 2: Ratio of simulated variance of the Horvitz-Thompson estimator for the two variables of interests. Two measures of spatial balance are displayed for each design. | $v_{sim}/v_{sim}^{CP}$ | Spread measures
---|---|---
| Foreigner percentage | Tertiary sector | $B$ | $I$
Only fixed-size
Equal
Proposed Method | 60.566 | 71.417 | 0.096 | -0.262
Local cube | 64.212 | 75.030 | 0.102 | -0.224
Local pivotal | 61.127 | 73.848 | 0.098 | -0.241
Max entropy | 100 | 100 | 0.383 | -0.003
Unequal
Proposed Method | 48.746 | 31.050 | 0.095 | -0.261
Local cube | 49.961 | 31.639 | 0.103 | -0.216
Local pivotal | 50.801 | 31.737 | 0.098 | -0.231
Max entropy | 100 | 100 | 0.390 | -0.003
Auxiliary variables
Equal
Proposed Method | 44.565 | 11.266 | 0.172 | -0.114
Local cube | 46.913 | 12.292 | 0.210 | -0.074
Local pivotal | 64.367 | 76.103 | 0.097 | -0.241
Max entropy | 100 | 100 | 0.382 | -0.003
Unequal
Proposed Method | 47.633 | 13.397 | 0.182 | -0.102
Local cube | 46.696 | 13.542 | 0.203 | -0.079
Local pivotal | 49.796 | 31.115 | 0.098 | -0.232
Max entropy | 100 | 100 | 0.390 | -0.003
Figure 2: Scatter plot of the auxiliary variables selected for the equal
inclusion probabilities design. Foreigners and Tertiary are variable of
interest. Figure 3: Percentage of foreigners and log of the number of
employment in tertiary sector for each municipality in Switzerland.

##  7 Conclusion

The increasing amount of data in our century will make sequential sampling a
powerful tool to handle big datasets. In some experiments in environmental
studies, it is even impossible to know the entire population and so sequential
sampling are mandatory. In this manuscript, we proposed a new sequential
method to select a balanced sample. This method respects equal and unequal
inclusion probabilities. We demonstrated through simulations that the proposed
method selects a more spread sample than the local cube method. Moreover, we
showed that the method significantly reduces the variance of the Horvitz-
Thompson estimator. To conclude, we can state that the method is an effective
improvement for sampling and a solid enhancement in this field of spatial
statistics.

## References

  * Chen et al., (1994) Chen, X.-H., Dempster, A. P., and Liu, J. S. (1994).  Weighted finite population sampling to maximize entropy.  Biometrika, 81:457–469.
  * Deville and Tillé, (2000) Deville, J.-C. and Tillé, Y. (2000).  Selection of several unequal probability samples from the same population.  Journal of Statistical Planning and Inference, 86:215–227.
  * Deville and Tillé, (2004) Deville, J.-C. and Tillé, Y. (2004).  Efficient balanced sampling: The cube method.  Biometrika, 91:893–912.
  * Gini and Galvani, (1929) Gini, C. and Galvani, L. (1929).  Di una applicazione del metodo rappresentativo al censimento italiano della popolazione (1. dicembre 1921).  Annali di Statistica, Series 6, 4:1–107.
  * Grafström and Lisic, (2019) Grafström, A. and Lisic, J. (2019).  BalancedSampling: Balanced and Spatially Balanced Sampling.  R package version 1.5.5.
  * Grafström and Lundström, (2013) Grafström, A. and Lundström, N. L. P. (2013).  Why well spread probability samples are balanced?  Open Journal of Statistics, 3(1):36–41.
  * Grafström et al., (2012) Grafström, A., Lundström, N. L. P., and Schelin, L. (2012).  Spatially balanced sampling through the pivotal method.  Biometrics, 68(2):514–520.
  * Grafström and Tillé, (2013) Grafström, A. and Tillé, Y. (2013).  Doubly balanced spatial sampling with spreading and restitution of auxiliary totals.  Environmetrics, 14(2):120–131.
  * Hájek, (1981) Hájek, J. (1981).  Sampling from a Finite Population.  Marcel Dekker, New York.
  * Jauslin et al., (2021) Jauslin, R., Eustache, E., Panahbehagh, B., and Tillé, Y. (2021).  StratifiedSampling: Different Methods for Stratified Sampling.  R Foundation for Statistical Computing, Vienna, Austria.  R package version 0.3.0.
  * Jauslin and Tillé, (2020) Jauslin, R. and Tillé, Y. (2020).  Spatial spread sampling using weakly associated vectors.  Journal of Agricultural, Biological and Environmental Statistics, 25(3):431–451.
  * Leuenberger et al., (2021) Leuenberger, M., Eustache, E., Jauslin, R., and Tillé, Y. (2021).  Balancing a sample almost perfectly.  Statistics and Probablity Letters, accepted.
  * Moran, (1950) Moran, P. A. P. (1950).  Notes on continuous stochastic phenomena.  Biometrika, 37(1/2):17–23.
  * Stevens Jr. and Olsen, (2003) Stevens Jr., D. L. and Olsen, A. R. (2003).  Variance estimation for spatially balanced samples of environmental resources.  Environmetrics, 14(6):593–610.
  * Stevens Jr. and Olsen, (2004) Stevens Jr., D. L. and Olsen, A. R. (2004).  Spatially balanced sampling of natural resources.  Journal of the American Statistical Association, 99(465):262–278.
  * Swiss Federal Statistical Office, (2021) Swiss Federal Statistical Office (2021).  Regionalporträts 2021: Kennzahlen aller Gemeinden.  Neuchâtel, Switzerland.
  * Tillé et al., (2018) Tillé, Y., Dickson, M. M., Espa, G., and Giuliani, D. (2018).  Measuring the spatial balance of a sample: A new measure based on the Moran’s ${I}$ index.  Spatial Statistics, 23:182–192.
  * Vallée et al., (2015) Vallée, A.-A., Ferland-Raymond, B., Rivest, L.-P., and Tillé, Y. (2015).  Incorporating spatial and operational constraints in the sampling designs for forest inventories.  Environmetrics, 26(8):557–570.