Noname manuscript No.
(will be inserted by the editor)

Fast Selection of Spatially Balanced Samples
F. Piersimoni · R. Benedetti

arXiv:1710.09116v1 [stat.ME] 25 Oct 2017

the date of receipt and acceptance should be inserted later

Abstract Sampling from very large spatial populations is challenging. The solutions suggested in recent literature on this subject often require that the randomly selected units are
well distributed across the study region by using complex algorithms that have the feature,
essential in a design–based framework, to respect the fixed first–order inclusion probabilities for every unit of the population. The size of the frame, N , often causes some problems
to these algorithms since, being based on the distance matrix between the units of the population, have at least a computational cost of order N 2 . In this paper we propose a draw–
by–draw algorithm that randomly selects a sample of size n in exactly n steps, updating at
each step the selection probability of not–selected units depending on their distance from the
units already selected in the previous steps. The performance of this solution is compared
with those of other methods derived from the spatially balanced sampling literature in terms
of their root mean squared error (RMSE) using the simple random sampling (SRS) without
replacement as benchmark. The fundamental interest is not only to evaluate the efficiency of
a such different procedure, but also to understand if similar results can be obtained even with
a notable reduction in the computational burden needed to obtain more efficient sampling
designs. Repeated sample selections on real and simulated populations support this perspective. An application to the Land Use and Land Cover Survey (LUCAS) 2012 data–set in an
Italian region is presented as a concrete and practical illustration of the capabilities of the
proposed sample selection method.
Keywords Spatial dependence · Big Data · Generalized Random Tessellation Stratified
design · Spatially correlated Poisson sampling · Local pivotal method · Product of the
within sample distance
.
F. Piersimoni
Istat, Directorate for Methodology and Statistical Process Design, Via Cesare Balbo 16, Rome, IT-00184,
Italy
E-mail: piersimo@istat.it
R. Benedetti
“G. d’Annunzio” University, Department of Economic Studies (DEc), Viale Pindaro 42, Pescara, IT-65127,
Italy
Tel.: +39-085-45083210, +39-085-45083229
Fax: +39-085-45083208
E-mail: benedett@unich.it

2

F. Piersimoni, R. Benedetti

Mathematics Subject Classification (2000) 62D05 · 62H11

1 Introduction
Surveys are routinely used to gather data for environmental and ecological research. The
units to be observed are often randomly selected from a population made up of geo-referenced
units, and the spatial distribution of these units is information that can be used in sample design. It represents a source of auxiliaries that can be helpful to design an effective random
selection strategy, which, by a proper use of this particular information and assuming that
the observed phenomenon is related with the spatial features of the population, could lead
to a remarkable efficiency gain of a design–based estimator of the unknown total of a target
variable.
In recent years, with the widespread use of the Global Positioning System (GPS) technology,
of remotely sensed data, of automatic address geocoding systems and of Geographical Information Systems (GIS) databases, spatial modeling and analysis of geo-referenced data-sets
have received more and more attention. The use of this particular type of data is becoming popular in several research areas including geological and environmental sciences [62],
ecology [51], cartography [37], public health [15] and so on. Some extended reviews of current research tendencies in spatial modeling and analysis are presented in [11, 19, 23, 25, 40,
42,50].
At the same time it is surprising and incomprehensible the reason why, except for some attractive references on the subject [6,44], at least the same interest was not devoted to basic
topics as spatial data collection and sample surveys selected from spatial populations, with
the consequence that none of these authoritative texts addresses, unless mentioned as a border issue, the problem of sampling from spatial data–sets.
To better understand the increase in the importance of spatial data, it is also appropriate
to realize that is more and more frequent the practice that National Statistical Offices georeference their sampling frames of physical or administrative bodies, used for social and
economic surveys, not only according to the codes of a geographical nomenclature, but
also adding information regarding the exact, or estimated, position of each record. With the
effect that such methodological problems are no longer felt only in local and occasional experiments, although of a high scientific level, but also in economic and social surveys that
provide basic statistical figures at national and regional level. Indeed area frame surveys [6]
often supported, or even replaced, the traditional methods based on list frames mainly due
to the absence of coverage errors and the low probability of non–responses that they guarantee. Moreover the latest fashion in official statistics, currently considered a keyword, is
the statistical use of administrative data that implies the need of estimating the coverage of
administrative lists, which is often carried out through area frame surveys. Spatial surveys
also present some drawbacks with regard to the reduced information that can be collected
in a survey involving only a direct observation and not a questionnaire to be filled. Within
this context, in [2] it is claimed that area frame surveys will play a very important role in the
future development of agricultural surveys in particular in developing countries.
In addition big spatial data–sets are very common in scientific problems, such as those involving remote sensing of the earth by satellites, climate-model output, small-area samples
from national surveys, and so forth. The term Big Data includes data characterized not only
by their large volume, but also by their variety and velocity, the organic way in which they
are created, and the new types of processes needed to analyze them and make inference
from them. The change in the nature of the new types of data, their availability, and the

Fast Selection of Spatially Balanced Samples

3

way in which they are collected and disseminated is fundamental. This change constitutes a
paradigm shift for survey research. There is great potential in Big Data, but there are some
fundamental challenges that have to be resolved before their full potential can be realized.
A large data–set entails not only computational but also theoretical problems as it is often
defined on a large spatial domain, so the spatial process of interest typically exhibits local
homogeneity over that domain. The availability of massive spatial data that is becoming
quite common emphasizes the need for developing new and computationally efficient methods tailored to handle such big data–sets. These problems have already been widely treated
from a modeling approach and there is a lot of literature on the subject [12, 13, 35, 57]. On
the contrary, the aspect of data collection, in particular the sample selection, has not yet been
fully addressed.
The situation becomes even more difficult to manage if we consider that these data have
complex structures. The use of points, lines, polygons or regular grids to represent spatial
units and their relationships is certainly more than an established practice, it is the way used
in geography to simplify reality. Each data type in this list is represented by different modes
and rules representing a major framework for the inherent complexity of the analysis of
spatial phenomena. This complexity is related to the stage of data collection, in that various situations can arise for recording and describing a given phenomenon under alternative
topologies.
Reflecting on the evidence that spatial sampling algorithms can be very slow when the data–
sets are made up of a large amount of data, our main purpose was to propose a fast method
for random selection of units that are well spread in every dimension that can therefore be
used in massive surveys without creating computational problems and that, by founding its
ground on an attempt to set the concepts used by different approaches to the problem, have
the flexibility deriving from a model–based (MB) approach and can be safely used even by
the most strong supporters of sampling from finite populations within a design–based (DB)
framework.
The plan of the paper is as follows. Some preliminaries and basic concepts are given in Section 2 together with some motivation to spread the sample over a spatial population. After
a review of the literature on spatially balanced samples, in Section 3 is introduced a fast
algorithm that, selecting well–spread samples, seeks to exploit the spatial characteristics of
the population. In Section 4 the performance of the proposed method is empirically compared with some reviewed competitors by a design–based simulation study. The efficiency
appraisal is evaluated in terms of the root mean squared error (RMSE) of the estimates by using the simple random sampling (SRS) as benchmark. Different estimation scenarios have
been assumed for varying sample sizes and with respect to real and artificial populations
characterized by different spatial distributions. Finally, in Section 5 the main results of the
paper are discussed and some addresses for further research are provided.

2 Basic concepts and preliminaries
The adoption of a spatial model exploiting the topological nature of geographical entities
can be formally tackled by introducing the following super–population model for spatial
data [11]: let i ∈ Rd be a generic data location in a d–dimensional Euclidean space within a
region U ⊆ Rd and Y i a random variable representing the potential datum for the attribute
or target variable at spatial location i. Now let i vary over index set U so as to generate the
random field {Y (i) : i ∈ U }. Within this model, space is featured in every respect by the set
U , including topology and coordinates system, and location in U indexed by i. Notice that

4

F. Piersimoni, R. Benedetti

the data generating process can lie also over a continuous spatial domain even if, to simplify
the problem, it is observed only in a selection, possibly made at random, of fixed points or
averaged over a selection of predefined polygons. Even if in this paper we will focus our attention on situations of this kind, which can be led back to usual finite population sampling,
it is appropriate to underline that in natural resources monitoring and estimation, however,
they cover an important part of all the possible sampling problems that arise in this field.
There is a huge list of phenomena that can be observed in any site of a linear object, such as
a river, or of a surface as it is for meteorological data. In these cases the resulting sample is
a set of points or polygons whose possible positions are not predefined but chosen from an
infinite set of possible sites.
A spatial finite population U = {1, . . . , i, . . . , N } of size N is thus recorded on a frame and
from the geographical position and topology of each unit i we can derive, according to some
distance definition, a matrix DU = {dij ; i = 1, . . . , N ; j = 1, . . . , N } that specifies how
far are all the pairs of units in the population.
The use of the matrix DU as a synthesis of the spatial information implies the hypothesis
that the dependence doesnot change with the position of the unit i and the direction, i.e.
that the random field Y i is homogeneous and isotropic [11], i.e. its distribution does not
change if we shift or rotate the space of the coordinates. DU is a very important tool to emphasize the importance to spread the sample over U, a property that can be related, through
a variogram [11], to the spatial dependence of Y i but also to some form of similarity between adjacent units as a spatial clustering or a spatial stratification. If the units are points,
for the definition of DU , we can simply resort to simple concepts of distance between sets
of coordinates, but if they are polygons we should use as a distance the notion of contiguity between areal units, or it would be better to use the order of contiguity, unless we want
to transform polygons into points by identifying them with their centroids. As far as linear
topology is concerned, a distance measurement can not ignore assumptions about the distances between the sets of points that make up the lines. Some synthetic indexes such as the
minimum and the maximum are always available while for the average, as it is necessary to
calculate the area between the two lines, serious computational problems

P could arise.
A traditional objective of most surveys is estimation of the
R totaltY = i∈U Y i when we
assume that U is a finite population that becomes tY = U Y i d i if U is not finite but is a
continuous surface over U . In order to estimate tY , a sample of units from U is selected by
identifying their labels i on the frame, and then measures yi of their corresponding values
are collected. Considering only without replacement samples, a convenient way to state this
sample selection process is to assume that, for each unit i on the frame, a random vector
S = {si ; i = 1, . . . , N } is generated that is equal to 1 if a unit is selected in the sample and

P
0 otherwise. If the sample size is fixed then n = N
i=1 si . The distribution P S , whose
support is S and in principle is under the complete control of the sampler, defines what is
generally referred to as the design of the sample survey. The first–order
inclusion probability
P
that the unit i will be included in the sample is denoted by πi = S3i P (S ), where the term
S 3 i means that the sum is extended over those samples that contain
P i, while the second
order inclusion probability for the units i and j is denoted as πij = S3{i,j} P (S ).
Within the DB estimation framework we consider P (S ) as the only source of randomness,
thus assuming that yi is not a random field and is not affected by any measurement error, we
can define the well-known, and widely usedP
in practical applications, Horvitz and Thompson (HT) [49] estimator for tY as t̂HT,Y = i∈S yi /πi , where the term i ∈ S here means
that the sum is extended over those units i that have been selected in the sample S ∈ S .
The equivalent of the HT estimator, and its properties, for continuous populations was ex-

Fast Selection of Spatially Balanced Samples

5

tensively treated by [14].
The alternative MB estimation framework consists in assuming that the actual
finite popu
lation yi is only one of the possible realizations of the random field Y i . The population
total can be decomposed as the sum of two components tY = tYS + tYS̄ , where YS is the set
of yi values observed in the sample while YS̄ are the not–observed values in the remaining
units of the population (S̄ ≡ U − S ). Once completed the data collection tYS is known, and
the estimation problem can be reduced to predict tYS̄ with t̂YS̄ as the sum of the predicted
values ŷi arising from a suitable model ξ fitted to the observed data [9, 58].
It follows that, depending on which estimation approach we adopt, we should evaluate the
expected value and the mean squared error (MSE) of the estimators following different crite


2 
ria, when dealing with a DB approach ES t̂HT,Y −tY and ES t̂HT,Y −tY
should be

 

2 
used while Eξ t̂HT,Y − tY |S and Eξ t̂HT,Y − tY |S are more appropriately applied
when working within a MB approach, where ES and Eξ denote respectively expectation
with regard to the sample S ∈ S and to the model ξ .
A vector of C covariates xi = {xi1 , . . . , xic , . . . , xiC } is usually available for each unit i
in the frame represented by the coordinates of the unit (if they are points or centroids of
polygons), the land use derived from a map, the elevation, remotely sensed data, administrative data and so on. Considering that we are dealing with spatial data the practice is that
the model ξ consists of two parts, a trend component that is used to relate the outcome Y i
with the auxiliaries xi and an autocorrelation component that takes advantage of the knowledge of DU to fit a variogram to the observed data. This model is known as kriging [11, 19]
and, considering its theoretical properties and its popularity in practical applications, it is
surely one the best candidate for ξ .
To select S the logic of the optimal design is conceived so that preferential sampling [18]
can be used allowing explicitly for the minimization of a criterion Φ linked with some summary statistics Q, usually an utility function, arising from model ξ [61]. More formally the
set of n units that constitute the sample S , and possibly a set of weights
 F associated to the
sample, are the result of an optimization problem such as max Q ξ [43, 45, 46]. This is a
S,F

combinatorial optimization problem that can be solved through some heuristics or by using
the well known Simulated Annealing algorithm that has shown to provide promising results
[4]. Alternatively a design criterion motivated from Bayesian learning was suggested and
estimated by MCMC [22] or, if estimation of Φ is complex as the likelihood is intractable,
approximate bayesian computation showed to be a feasible solution [31, 32].
In a MB perspective, the concern is necessarily in finding the sample configuration that
is the best representative of U maximizing an objective function defined over the whole
set S of possible samples which, in a spatial context, will surely depend on the loss of
information due to spatial autocorrelation. The resulting optimal sample is selected with
certainty and is of course not advised, if we assume the randomization hypothesis that is the
background for DB inference [48,56,60]. Moreover, as clearly stated by [12]: “A modelbased approach seems therefore necessary under preferential sampling, such as assuming a
spatial-statistical model for s(·). However, this does not mean that the basic design-based
notions of randomisation, stratification, and clustering cannot be used in a preferential–
sampling approach, since they are all useful tools that lead to a better representation of a
heterogeneous population. In particular, what happens when the model . . . does not adequately describe the spatial variability in y (·) and s(·)? The optimality . . . , and the validity
of any consequent inference depends critically on the appropriateness of this model”.
The answer to these questions is, of course, a highly debated topic not only for theoretical
but also practical reasons, starting with the evidence that if we do not assume that a model

6

F. Piersimoni, R. Benedetti

holds for the measurement errors, the potential observations over each unit of the population, being of a deterministic nature, cannot be considered as dependent. Thus it is clear that
sampling schemes to be reasonably adopted for spatial units in a DB context need a suitable
alternative to the concept of spatial dependence.
A useful property of random samples, that of being spatial balanced, has been introduced
for this purpose by [53]. To define a spatial balance index (SBI) the use of Voronoi polygons
is required for any sample S . Such a polygon, for the sample unit i, includes all population
units closer to i than to any other sample unit j . If we let νi be the sum of the πv of all units
v in the i-th Voronoi polygon, for any sample unit we have E (νi ) = 1. Thus the SBI:
P
SBI (S ) = V ar (νi ) =

i∈S

(νi − 1)2
n

,

(1)

can be used as a measure of difference from the state of perfect spatial balance (SBI(S)=0).
A sample S is considered spatially balanced when SBI(S) ≤ ωsb as every νi should be close
to 1, where ωsb is an acceptable upper bound. Note that the role played by SBI is exactly
the same as Φ, it is nothing more than a criterion for assessing the suitability of S but,
depending on the πi s, its use is more appropriate in a design–based framework. Anyway, its
minimization would again lead to a purposive sample that would not solve the problem of
randomization of S .
A way to recover the use of autocorrelation in a DB sampling strategy, or at least to restore
its acceptability, is to consider that in survey design an upper bound to the sampling error is
usually imposed to the auxiliary variables rather than on the target variable. This choice is
usually motivated by the practitioners on the basis of the hypothesis that a survey that places
its bases on specific levels of precision for a set of auxiliary variables will approximately
yield the same sampling errors also for the target variable. However, in realistic cases of
practical interest, this assumption is unlikely
to occur and appreciable differences are often

measured among the covariates and Y i. In such situations the suggested design could be
incorrect as using xi as a proxy for Y i it could underestimate the sample size needed to
reach a predetermined level of precision.

An appealing and widely used alternative is to assume a model that links Y i to xi . The
assumption underlying this solution is that it is possible to estimate the unknown parameters
of such models from data collected in previous surveys. This concept of anticipated variance
(AV - or anticipated MSE if the estimator is biased) was introduced by Isaki and Fuller [34]
and is very useful to overcome the hitch of not being able to introduce, at least in the design
of the sample, a stochastic model that can take into account the dependence of spatial units.
It is defined as the average of the DB variance of t̂HT,Y , the HT estimator of the total of
Y , under a stochastic model ξ . To predict the Y i s we assume that, for each unit i of U , a
linear model ξ holds, given the known auxiliaries xi s:


Y i = xti β + i



Eξ (i ) = 0
V
arξ (i ) = σi2



Covξ (i j ) = σi σj ρij

,

(2)

where Eξ , V arξ , and Covξ denote respectively expectation, variance and covariance with
respect to the model ξ , β is a vector of regression coefficients, i is a random variable with
variance σi2 and ρij is its autocorrelation coefficient. The AV of the unbiased HT estimator

Fast Selection of Spatially Balanced Samples

7

of the total of Y , under model (2), is:

!T 
X xi X
XX

πij − πi πj
AV t̂HT,Y − t = ES 
−
xi
β +
.
σi σj ρij
i∈S

πi

i∈U

i∈U j∈U

πij

(3)

The two components of (3) represent respectively the error implicitly introduced in the estimate of the totals of the auxiliary variables and the existing dependence of the population
units. It is clear that uncertainty on estimates can be reduced constraining the units selected
to respect the average value of the population’s covariates and, assuming that autocorrelation
coefficient decreases as the distance dij between the selected units i and j increases, selecting units as far apart as possible. The logic derived from the AV criterion may, however,
sound like an excuse to justify the forced introduction of a model ξ within a DB framework.
In this regard, it is interesting to note that the authors of this proposal justify the designs derived from it within a MB and not DB framework [29]: “Even though we justify the method
by using a superpopulation model, the inference is based on the sampling design. The HT
estimator will be efficient if the population is close to a realization from the model, but the
estimator maintains desirable properties like design unbiasedness and design consistency
even if the model is not properly specified”.
Another motivation for spatially balanced samples, perhaps less forced and even less arbitrary, has been recently proposed [5] on the basis of the so called decomposition lemma, a
useful though not so well-known result in sampling theory, which states that [39, p. 87]:
¯S +
σy̆2 = VS y̆


n−1
ES s2y̆,S ,
n





(4)

where y̆ is a vector of the expanded-values of the target variable y whose generic element is
¯S
y̆i = πi−1 yi , σy̆2 is the constant and unknown population variance of the variable y̆, VS y̆
¯S = 1 P
y̆
is the variance between samples of the HT estimator of the mean y̆
i∈S i and
n

P
2
1
¯
ES s2y̆,S is the expectation of s2y̆,S = (n−
(y̆
−
y̆
)
,
i.e.
the
within
sample
variS
i∈S i
1)
ance according to the design P (S ) (for details see [39, ch.3]).
It can be seen from (4) that a gain in the efficiency of the HT estimator can be realized either
by setting the πi s in such a way that y̆ is approximately constant [18, p. 53] or by defining
a design P (S ) that increases the expected within sample variance or both. Provided that we
are dealing with samples S with fixed and known πi s, the intuitive explanation for this is
that, being constant the unknown variance of the population, the only way to reduce the first
term in the right side of (4) is to increase the second term. Thus, we should set the probability P (S ) to select a sample S proportional, or even more than proportional, to s2y̆,S .
Unfortunately, this proposal risks to remain purely theoretical as this parameter is unknown
being relative to the unobserved target variable y. In the spatial interpolation literature [11,
19], it is often adopted the assumption that the distance DU is highly related to the variance
of a variable observed on a set of geo–referenced units. The variogram (or semi variogram),
and its shape, is a classical tool to choose how and to what extent the variance of y̆ is a
function of the distance between the statistical units. Thus, when dealing with spatially distributed populations a promising candidate to play the role of the proxy for s2y̆,S is DU . The
intuitive requirement that a sample should be well–spread over a study region is explainable
if and only if there are reasons to assume that s2y̆,S is a monotone increasing function of
DU . This will surely happen when y̆ has a linear or monotone spatial trend or when there is
spatial dependence. It follows that the distance
 matrix DS between sample units, if appropriately synthesized through an index M DS , could be an additional criterion to play the

8

F. Piersimoni, R. Benedetti

same role as SBI and Φ and is probably much simpler and faster to evaluate.
Finally, noting that the Sen-Yates-Grundy (SYG) formulation of the estimator of the sample
variance:



1 X X πij − πi πj
2
V̂SY G t̂y = −
(5)
(y̆i − y̆j ) ,
2
πij
i∈S j∈S

is very similar to the definition of the variogram, another interesting motivation to the practice to spread the sample can be derived from the assumption that, increasing the distance
2
dij between two units i and j , the difference (y̆i − y̆j ) between the values of the survey
variable is expected to increase. As a result an efficient design, to weigh less the expected
most relevant differences (y̆i − y̆j )2 , should use a set of πij strictly and positively related
with dij .

3 Algorithms for random selection of spatially balanced samples
The need to avoid as much as possible to include in the same sample contiguous or too close
spatial units is very much felt in real periodic surveys such as those used to estimate land use
and land cover [6, ch. 2] or to implement a forestry inventory [30, 41]. The strategies adopted
in planning these surveys invariably invoke simple and consolidated methods of sampling
theory that have the drawback of not being specifically implemented to solve this problem.
These solutions often involve systematic selection with a predetermined step and a random
starting point, partitioning the area in exactly n contiguous strata from which randomly pick
up only one unit per stratum and multi-stage samples with primary units represented by
spatial aggregates. Each of these proposals has proved not to be appropriate to capture the
potential efficiencies deriving from the existence of spatial trends or dependence. The main
reason for these difficulties lies in the evidence that, if the population does not lie on a regular grilling, it is almost impossible to define a sampling step in the absence of a unit order
and it is highly subjective to set up a partition of the region in n parts. The choices made to
overcome these shortcomings seriously affect the results.
One of the first attempts to formalize the problem, increasing the amount of information
collected by avoiding the selection of pairs of contiguous units, was done in the pioneering
work [33] that, however, explicitly requires the introduction of an exogenous ordering of the
population units. The line drawn from this preliminary work has remarkably influenced the
approaches adopted in the following studies on this topic. Many algorithms, even recently
proposed in the literature [53], are primarily based on the search for a one-dimensional sorting of multidimensional units by studying the best mapping of Rd in R, while trying to
preserve some multidimensional features of the population, and then use this induced ordering to systematically select the sample units. The fundamental principle is to extend the
use of systematic sampling to two or more dimensions, even when the population is not a
regular grid. Developed by the Environmental Protection Agency (EPA) and widely used in
most of its surveys, the Generalized Random Tessellation Stratified (GRTS) design [53] is
mainly based on such a logic. A multilevel grid hierarchy is used to index the units through
a tree structure.
The two parts of the AV (3) suggested instead the use of two different classes of methods, the
first derived from the restriction of the sample space considered acceptable and the second
from the assumption that dependence between units decreases as their distance increases.
Despite the assonance of the two names, the notion of spatially balanced samples is quite far
from that of balanced samples used to denote those samples that respect the known totals of

Fast Selection of Spatially Balanced Samples

9

a set of auxiliary variables, their selection is made at random on a restricted support S̄ ⊂ S of
all the possible samples [54,55,56]. Thus, for given πi s the set of balanced samples are defined as all those samples S for which it occurs |t̂HT,xc − tHT,xc | ≤ ωb , ∀ c = {1, . . . , C},
where ωb is an acceptable upper bound. These restrictions represent the intuitive requirement that the sample estimates of the total of a covariate should be as close as possible to
the known total of the population. In a spatial context, this constraint could be applied by
imposing that, for any S , the first p moments of each coordinate
 should coincide with the
first p moments of the population, implicitly assuming that Y i follows a polynomial spatial trend of order p. This logic was subsequently extended to approximate any nonlinear
trends through penalized splines with particular reference to the space [8].
The CUBE algorithm [10,17], proposed as a general (not spatial) solution, is the tool used
to achieve these samples. In recent studies [5,6], these sampling designs have demonstrated
that they can effectively exploit the presence of linear and nonlinear
 trends but are often
unable to capture the presence of dependence or homogeneity in Y i .
The second class is instead more articulated and consisting of algorithms of different nature.
A selection strategy conceived with the clear goal of minimizing (1) should use DU as the
only known information available on the spatial distribution of the sample units. Its use implies the adoption of the already mentioned intuitive criterion that units that are close should
seldom appear simultaneously in the sample. Recalling (5), this request can be claimed as
reasonable if we expect that, increasing dij , the difference (y̆i − y̆j )2 always increases. In
such a situation, it is clear that the variance of the HT estimator will necessarily decrease if
we set high joint inclusion probabilities to couples with very different values of the target
variable as they are far each other. Following an approach based on distances, inspired by
purely MB assumptions on the dependence of the stochastic process generating the data,
[1] suggested a sampling strategy: the dependent areal units sequential technique (DUST).
Starting with a unit selected at random, say i, at every step t < n, the selection probabilities
are updated according to a multiplicative rule depending by a tuning parameter useful to
control the distribution of the sample over the study region. This algorithm, or at least the
design that it implies, can be easily interpreted and analyzed in a DB perspective in particular referring to a careful estimation and analysis of its first and second order inclusion
probabilities.
Another solution, based on a classic list sequential algorithm [54], was suggested by [26].
Introduced as a variant of the correlated Poisson sampling, in the SCPS (Spatially Correlated Poisson Sampling) for each unit, at every step t, it updates the inclusion probabilities
according to a rule in such a way that the required inclusion probabilities are respected. The
suggested maximal weights criterion, used to update the inclusion probabilities at each step,
provides as much weight as possible to the closest unit, then to the second closest unit and so
on. The resulting sample is then obtained in n ≤ t ≤ N steps depending on the randomness
according to which the last unit is added.
A procedure to select samples with fixed πi and correlated πij was derived by [28] as an extension of the pivotal method initially introduced to select π ps samples [16]. It is essentially
based on an updating rule of the probabilities that at each step should locally keep the sum
of the updated probabilities as constant as possible and differ from each other in a way not
to choose the two nearby units. This method is referred to as the local pivotal method (LPM)
and obtain a sample in, again, n ≤ t ≤ N steps. This solution has been also integrated with
the balancing property, typical of the CUBE algorithm, to simultaneously minimize both the
two members of the AV (3) [29]. A comprehensive review of the main spatially balanced
samples selection methods can be found in [7].
According to (4) the attention should be moved to the definition of a design with probability

10

F. Piersimoni, R. Benedetti


proportional to some synthetic index M · of the within sample
distance matrix when it
P
is observed within each possible sample: P (S ) = M DS / S∈S M DS [5]. This innovative proposal is quite dangerous because extracting realizations from the complete P (S )
could make it very difficult to control or specify the πi s and the πij s whose knowledge is
required by the HT estimator to be able to yield total and variance estimates. Despite this
initial difficulty it may open a very interesting window on how to use indices, such as Φ or
SBI, within a DB framework. Select the whole sample with probability proportional to the
index instead of maximizing or minimizing the index itself could remove all the doubts that
sometimes arise from choosing the units with certainty. Interpreted in this way, this proposal
could thus represent a bridge between the MB and DB frameworks to setting up a spatial
sampling strategy.

Among the possible summary indexes
· of the distance matrix, the most promising reQMQ
sults were provided by M (DS ) =
dγij introducing a design proportional to the
i∈S j6=i;j∈S

products of the within sample distance matrix (PWD) that depends on γ , a tuning parameter
that can be used without limits to increase or decrease the effects of the distance matrix on
the spread of the sample over the study region. In [5] it is suggested to start with a SRS
without replacement and repeatedly and randomly exchanging a unit included in the sample
with a unit not included in the sample with probability equal to the exponential of the ratio
between the two indexes before and after the exchange. This is an MCMC iterative procedure that for a suitable choice of the number of iterations will generate a random outcome
from the multivariate distribution P (S ).
These procedures are quick and efficient and thus have a practical applicability in real spatial
surveys if the size N of U is not prohibitive. The main reason is that they strictly depend on
N in the number of attempts or steps needed to select a sample. According to the definition
[54, p. 35] “A sampling design of fixed sample size n is said to be draw–by–draw if, at each
one of the n steps of the procedure, a unit is definitively selected in the sample”, a possible
way to reduce the computational burden is to look for a draw–by–draw alternative to the
listed methods. With regard to the PWD, deriving its probabilities from the product of the
distance between units, a natural candidate could arise from the idea of iteratively updating,
through a product, the πi s at each step.
The algorithm suggested in this paper, starts by randomly selecting a unit i with equal probability. Then, at every step t ≤ n, the algorithm updates the selection probabilities πjt s of
every other unit j of the population according to the rule (see Figure 1):
πjt = P

πjt−1 d¯ij

j∈U

πjt−1 d¯ij

∀j ∈ U,

(6)


where d¯ij = Ψ dγij is an appropriate transformation applied toQDU in order to standardize it so that it will have known and fixed products by row i6=j,i∈U dij and column
Q
t
¯
i6=j,j∈U dij . Notice that, as dii = 0 by the definition of distance, πi will be necessarily equal to 0 for all units i already selected in the sample in previous steps, thus preventing
random selection with replacement. Criterion (6) can be basically considered as a heuristic
method (HPWD) to generate samples approximately with the same probabilities of the PWD
but with a much smaller number of steps. From a practical point of view, the idea behind (6)
is very similar to the DUST method already suggested in the literature [1] although with two
substantive differences: (6) is motivated by solid theoretical bases while the DUST did not
find any justifications either according to a MB or DB logic and in the HPWD it is possible
to control the πi s while with DUST this is impossible with the result that the HT estimator

Fast Selection of Spatially Balanced Samples

11

Fig. 1 First order inclusion probabilities πjt of the HPWD of a population on a 100 × 100 regolar grid
after the selection of the first t = 1, . . . , 6 units with the suggested algorithm (darker is lower and lighter is
higher).

systematically yields irreparably biased estimates.
In [5] it is advised that the πi s and the πij s of the PWD follow the rules:

π̂i = k1 

N
Y

k2
dij 

(7)

,

j =1


π̂ij = k3 

N
Y

i=1

dij ×

N
Y

k4
dij 

(dij )k5 ,

(8)

j =1

where k1 , k2 , k3 , k4 and k5 are parameters mainly depending on the sample and population
sizes. Notice that, to ensure the symmetry of the πij s, the rows and columns marginal products of the distance matrix have the same parameter.
From (7) it is clear that if we want to fix the πi s, it is enough to standardize the dij s using the
d¯ij s specifying an appropriate standardizing function Ψ . For this role [5] suggest to iteratively constrain to known totals, the rows (or columns) sums of the logarithmic transformed
matrix. When the known totals are all equal to a constant, this is known to be a very simple
and accurate method to scale a symmetric matrix to a doubly stochastic matrix [36].
To verify if (7) and (8) can be used as working rules to specify the πi s we can generate as
many independent replicates from HPWD as needed and the πi s and the πij s may be estimated by the frequency with which a unit or a pair of units are selected. These π̂i s and π̂ij s
can both be adopted in the estimation process instead of their theoretical counterparts [20,
21] and also modeled to verify the fit of (8).

12

F. Piersimoni, R. Benedetti

Table 1 Relative efficiency of the estimated first order inclusion probabilities CV (π̂i ) (9) of the HPWD
estimated in 100,000 replicated samples in the 5 × 5 and 10 × 10 regular grid populations for different
sample sizes and γ.
5×5
10 × 10
γ n = 5 n = 10 n = 5 n = 10 n = 20 n = 40
1 1.852
0.653 2.864
2.688
1.477
0.453
5 8.543
2.457 13.373
8.082
3.659
1.469
10 12.722
5.297 21.358 13.131
6.080
2.583

To this purpose we selected 100,000 replicated samples of size n = {5, 10} from a 5 × 5
regular grid and of size n = {5, 10, 20, 40} from a 10 × 10 regular grid by using the HPWD
design with 3 different values of the parameter γ ={1, 5, 10}.
As far as the empirical frequencies of units in selected samples are concerned, in Table (1)
are shown the values of the index:

CV (π̂i ) =

v
uP


n 2
u
π̂
−
t
i
i∈U
N
N

n

N

× 100.

(9)

From the results in Table (1) we can observe that the πi s are approximately constant once
the row (and column) products of the d¯ij s are fixed as constant with an error that increases
with γ and decreases with the increase of n. The first effect is probably due to the fact that
increasing γ necessarily increases the concentration of the d¯ij s in few very high values and
many low values, which undoubtedly reduces the accuracy of the matrix standardization and
of the logarithmic transformation. While the second aspect occurs as increasing the number
of random selections the empirical π̂i s are closer to the theoretical πi s.
With regard to the π̂ij s the evidence from Figure (2) and Table (2) is that (8) fits enough
well supporting the assumption of such a relationship between the distance and the estimated
second order probabilities. However, some trends can be noticed not only in R2 but also in
the intercept (k̂3 ) and in the slope (k̂5 ) estimated through ordinary least squares (OLS) on
a logarithmic transformation of (8). As expected increasing γ the effects of the distance on
the π̂ij s are more sensible leading to a more scattered spatial distribution of the sampling
units. An evidence that is mitigated by the increase of n as in this case the HPWD has less
space to better spread the units.
Regarding the evaluation of the πij s we should also consider that their use is expected
only in HT variance estimation and that a spatial design producing well spread samples
necessarily will have very small joint inclusion probabilities for nearby units. Even though
HPWD has the advantage to produce πij s that can be shown to be always strictly positive,
although it is still unbiased, the HT variance estimator could become extremely unstable if
there is some very small πij . Thus, using HPWD, the HT estimation of variance is always
possible, and this is a clear advantage, but it cannot always be recommended as a good
solution. In such circumstances it could be better the use of alternatives variance estimators
as a MB approach to the problem [3] or the local variance estimator suggested by [52] that
represents a general solution as it does not need the πij s.

Fast Selection of Spatially Balanced Samples

13

Fig. 2 Second order inclusion probabilities π̂ij of the HPWD estimated in 100,000 replicated samples in the
5 × 5 and 10 × 10 regular grid populations for different sample sizes and γ = 1.
Table 2 Regression parameters of the model (8) for the π̂ij of the HPWD estimated in 100,000 replicated
samples in the 5 × 5 and 10 × 10 regular grid populations for different sample sizes and γ.
population grid
5×5
5×5
5×5
5×5
5×5
5×5
10 × 10
10 × 10
10 × 10
10 × 10
10 × 10
10 × 10
10 × 10
10 × 10
10 × 10
10 × 10
10 × 10
10 × 10
∗∗∗

γ
1
5
10
1
5
10
1
5
10
1
5
10
1
5
10
1
5
10

n

 
log k̂3

k̂5

R2

5
5
5
10
10
10
5
5
5
10
10
10
20
20
20
40
40
40

−3.377(0.006)∗∗∗

0.718(0.012)∗∗∗

0.943
0.846
0.806
0.872
0.709
0.660
0.917
0.830
0.763
0.862
0.774
0.599
0.782
0.591
0.540
0.722
0.406
0.337

−3.025(0.045)∗∗∗ 3.022(0.085)∗∗∗
−2.668(0.096)∗∗∗ 5.563(0.180)∗∗∗
−1.842(0.005)∗∗∗ 0.388(0.010)∗∗∗
−1.606(0.027)∗∗∗ 1.205(0.051)∗∗∗
−1.426(0.048)∗∗∗ 1.918(0.091)∗∗∗
−6.149(0.003)∗∗∗ 0.831(0.004)∗∗∗
−4.951(0.052)∗∗∗ 5.799(0.073)∗∗∗
−4.116(0.074)∗∗∗ 11.453(0.105)∗∗∗
−4.587(0.003)∗∗∗ 0.598(0.004)∗∗∗
−4.066(0.014)∗∗∗ 2.235(0.020)∗∗∗
−3.039(0.052)∗∗∗ 5.485(0.074)∗∗∗
−3.157(0.002)∗∗∗ 0.361(0.003)∗∗∗
−2.916(0.010)∗∗∗ 0.988(0.014)∗∗∗
−2.714(0.017)∗∗∗ 1.607(0.024)∗∗∗
−1.789(0.001)∗∗∗ 0.165(0.002)∗∗∗
−1.725(0.004)∗∗∗ 0.316(0.006)∗∗∗
−1.685(0.007)∗∗∗ 0.430(0.010)∗∗∗

p < 0.001, ∗∗ p < 0.01, ∗ p < 0.05

Table 3 Relative efficiency of the estimated first order inclusion probabilities CV (π̂i ) (9) of the HPWD in
10,000 replicated samples in the LUCAS - Emilia Romagna population, Meuse river, clustered and sparse
simulated populations for different sample sizes and γ.

LUCAS - Emilia Romagna
Meuse river
Sim. Clustered Population
Sim. Sparse Population
γ n = 100 n = 300 n = 600 n = 10 n = 50 n = 50 n = 100 n = 200 n = 50 n = 100 n = 200
1
7.658
4.480
3.001 5.695 2.271 10.542
7.616
5.788 5.766
5.396
5.153
5 11.183
7.505
5.519 18.393 12.329 28.539 28.146 26.747 21.514 26.076 26.237
10 16.264 11.682
8.864 32.002 20.952 43.558 47.529 40.981 37.762 44.454 41.986

14

F. Piersimoni, R. Benedetti



RM SE
Table 4 Relative efficiency of the sample mean RM
and average SBI for each design estimated
SESRS
in 10,000 replicated samples in the clustered and sparse populations for different sample sizes, trend and
autocorrelation.

Design
GRTS
LPM
SCPS
PWD1
PWD5
PWD10
HPWD1
HPWD5
HPWD10
GRTS
LPM
SCPS
PWD1
PWD5
PWD10
HPWD1
HPWD5
HPWD10
GRTS
LPM
SCPS
PWD1
PWD5
PWD10
HPWD1
HPWD5
HPWD10

n
50
50
50
50
50
50
50
50
50
100
100
100
100
100
100
100
100
100
200
200
200
200
200
200
200
200
200

Clustered Population
Sparse Population
No Trend
Linear Trend
No Trend
Linear Trend
Autocorrelation Autocorrelation Average Autocorrelation Autocorrelation Average
Medium High Medium High SBI Medium High Medium High SBI
0.963 0.612 0.519 0.349 0.308
0.990 0.777 0.549 0.520 0.132
0.935 0.591 0.498 0.322 0.265
0.991 0.698 0.522 0.449 0.094
0.935 0.592 0.498 0.336 0.279
0.999 0.687 0.515 0.437 0.093
0.946 0.624 0.515 0.359 0.272
1.002 0.774 0.542 0.509 0.136
0.957 0.491 0.460 0.246 0.135
1.004 0.619 0.503 0.402 0.065
0.975 0.466 0.443 0.221 0.104
0.991 0.566 0.487 0.385 0.050
0.951 0.611 0.507 0.339 0.269
1.003 0.760 0.524 0.478 0.123
0.935 0.504 0.470 0.262 0.166
1.005 0.633 0.502 0.393 0.067
0.955 0.466 0.471 0.247 0.147
1.002 0.597 0.494 0.382 0.058
0.924 0.557 0.481 0.299 0.230
0.959 0.669 0.514 0.425 0.142
0.904 0.521 0.462 0.280 0.191
0.965 0.629 0.501 0.385 0.105
0.908 0.517 0.474 0.281 0.203
0.968 0.630 0.508 0.387 0.104
0.919 0.565 0.480 0.315 0.213
0.981 0.704 0.514 0.446 0.149
0.907 0.443 0.430 0.225 0.107
0.941 0.538 0.477 0.354 0.078
0.907 0.416 0.421 0.214 0.093
0.935 0.515 0.474 0.342 0.068
0.925 0.542 0.467 0.293 0.205
0.973 0.674 0.510 0.417 0.136
0.888 0.444 0.440 0.236 0.121
0.957 0.549 0.474 0.346 0.080
0.894 0.417 0.432 0.222 0.103
0.948 0.510 0.471 0.339 0.070
0.900 0.505 0.454 0.253 0.199
0.947 0.613 0.491 0.387 0.169
0.874 0.466 0.429 0.233 0.164
0.938 0.551 0.478 0.336 0.132
0.877 0.459 0.425 0.234 0.172
0.941 0.550 0.473 0.333 0.135
0.895 0.512 0.448 0.263 0.202
0.939 0.618 0.483 0.389 0.179
0.829 0.395 0.400 0.197 0.123
0.889 0.473 0.450 0.285 0.109
0.826 0.389 0.401 0.191 0.118
0.887 0.469 0.452 0.280 0.104
0.890 0.483 0.436 0.249 0.191
0.947 0.589 0.486 0.362 0.166
0.836 0.400 0.399 0.196 0.125
0.902 0.470 0.456 0.284 0.109
0.810 0.386 0.385 0.182 0.112
0.864 0.443 0.452 0.268 0.097



RM SE
Table 5 Relative efficiency of the sample mean RM
and average SBI for each design estimated
SESRS
in 10,000 replicated samples in the Meuse river population for two different sample sizes.
Sample Size n = 10
Sample Size n = 50
Design
Cadmium Copper Lead Zinc Av. SBI Cadmium Copper Lead Zinc Av. SBI
GRTS
0.897
0.887 0.884 0.908 0.173
0.839
0.809 0.806 0.838 0.176
LPM
0.884
0.869 0.862 0.901 0.129
0.762
0.724 0.736 0.757 0.133
SCPS
0.883
0.864 0.853 0.893 0.120
0.755
0.708 0.705 0.735 0.138
PWD1
0.907
0.895 0.905 0.915 0.178
0.801
0.770 0.774 0.797 0.185
PWD5
0.770
0.782 0.773 0.799 0.086
0.688
0.631 0.647 0.672 0.119
PWD10
0.704
0.747 0.727 0.747 0.059
0.664
0.602 0.619 0.646 0.110
HPWD1
0.902
0.888 0.919 0.915 0.181
0.798
0.760 0.775 0.805 0.189
HPWD5
0.767
0.782 0.784 0.797 0.085
0.678
0.621 0.648 0.671 0.125
0.722
0.754 0.735 0.754 0.058
0.671
0.609 0.626 0.649 0.112
HPWD10

Fast Selection of Spatially Balanced Samples

15

Fig. 3 Spatial distribution of the simulated populations: (a) clustered and (c) sparse. The variograms are
reported for Medium and High autocorrelation without trend (b) and for No Trend and Linear Trend without
autocorrelation (d).

4 Sampling Designs Comparison on Artificial and Real Populations
In this Section the performance of the HPWD is empirically compared with respect to alternative spatially balanced designs via 10,000 sample replications, which have been carried
out by using the free software environment for statistical computings R [47]. In particular,
we used the following R packages: BalancedSampling [27], and spsurvey [38]. For
further details on the R codes that can be employed to randomly select spatially balanced
samples, see [6, ch. 7].
As possible alternatives to the HPWD, we considered the GRTS [53], the SCPS [26], the
LPM [28] and the PWD [5]. In addition, for PWD and HPWD, three values for the parameter γ ∈ {1, 5, 10} were tried, indicated as suffix after the acronym of the design. To be
comparable each other, all these alternative designs have been used setting the first order
inclusion probabilities constant and equal to n/N .
The comparison between different designs has been performed by using the (RMSE) of the
HT estimates as relative to the RMSE obtained when using a SRS design that is used, thus,
as a scale factor to remove the known effects of the size of the population N and of the
sample n on the sampling errors. It is worth noticing that in every simulation performed, as
the HT estimator is unbiased, the RMSEs were always very close to the standard error of

16

F. Piersimoni, R. Benedetti

Fig. 4 Spatial distribution of the Meuse river data (a). The fitted exponential variogram models are reported
for the variables Cadmium (b), Copper (c) and Zinc (d).

each design as the bias can be considered negligible.
The experiments carried out concern both real and artificial populations. The latter are two
frames of size N =1,000 generated through point processes with two different levels of clustering of the units to control the distribution
of the coordinates and with different spatial

features of the response variable Y i . The coordinates {x1 , x2 } are generated in the unit
square [0, 1]2 according to a Neyman-Scott process with intensity of the cluster centers
equal to 10 [59] with 100 expected units per cluster. Two different scale parameters for
cluster kernel were used {0.005, 0.03}, representing respectively an highly clustered and
a sparse
population of spatial units (see Figure 3). For each frame four possible outcomes

Y i have been generated according to a Gaussian Markov Random Field with or withouta
spatial linear trend x1 + x2 + η , that explain approximately the 80% of the variance of Y i .
Spatial dependence of the errors η is modeled by an exponential variogram with two intensities of the parameter {0.01, 0.1}, representing respectively medium and high dependence
between units. To avoid the possible effects due to different variability, each population was
finally standardized to the same mean µY = 5 and standard deviation σY = 1.
To verify if n has any effect on the efficiency, from each population, samples of different

Fast Selection of Spatially Balanced Samples

17

Fig. 5 Land Cover LUCAS data on a 2km × 2km population grid in the Italian region Emilia Romagna
overlaid on the NUTS2 and NUTS3 boundaries (wider and narrower lines respectively).

size n ∈ {10, 50, 200} have been selected. One of the two real populations used is the well
known case study, largely debated and analyzed in the field of geostatistics, the Meuse River
data–set available in the package gstat. It is a set of 155 samples of top soil heavy metal
concentrations (ppm) used in our experiments as a population, along with a number of soil
and landscape variables. The samples were collected in a flood plain of the river Meuse,
near the Stein village (The Netherlands). Heavy metal concentrations are bulk sampled from
an area of approximately 15m × 15m. In addition to the topographical map coordinates the
top soil concentrations of four metals have been used: cadmium, copper, lead and zinc. On
these variables there is an high spatial autocorrelation as it is apparent from the variograms
(Fig. 4). In this population the experiment consists of selecting samples of size {10, 50}.
Land is very important for most biological and human activities, it is the main economic
resource for agriculture, forestry, industries, and transport. The information collected on the
land deals mainly with two interconnected concepts: land cover that refers to the biophysical
coverage of land (e.g. crops, grass, broad-leaved forest, or build-up area) and land use that
specifies the socio-economic use of land (e.g. agriculture, forestry, recreation or residential
use). It is for this reason that land use and land cover data collected by the EU survey LUCAS has been chosen as the second real population.
The Land Use/Cover Area frame Survey (LUCAS) [24] is a task by EUROSTAT that was
initially designed to deliver, on a yearly basis, European crop estimates and then became a
more general agro-environmental survey.
Our population has been built considering 2012 as reference year and the region Emilia–
Romagna in Italy as the area under investigation (see Figure 5). The regular grid 2km× 2km
of points constitutes our frame population. In the LUCAS experiment three different sample
sizes have been adopted {100, 300, 600} on data aggregated in 10 land cover classes.
The first check that must necessarily be done before performing the estimates is related to the

18

F. Piersimoni, R. Benedetti



RM SE
Table 6 Relative efficiency of the Land Cover codes area estimates RM
and average SBI for
SESRS
each design estimated in 10,000 replicated samples in the Emilia-Romagna LUCAS population for different
sample sizes.

Design
n
GRTS
100
LPM
100
SCPS
100
PWD1
100
PWD5
100
PWD10 100
HPWD1 100
HPWD5 100
HPWD10 100
GRTS
300
LPM
300
SCPS
300
300
PWD1
PWD5
300
PWD10 300
HPWD1 300
HPWD5 300
HPWD10 300
GRTS
600
LPM
600
SCPS
600
PWD1
600
PWD5
600
PWD10 600
HPWD1 600
HPWD5 600
HPWD10 600
Area (hec.)

Artif.
Land
0.981
0.956
0.958
0.965
0.956
0.951
0.959
0.947
0.948
0.935
0.908
0.922
0.928
0.900
0.892
0.945
0.904
0.895
0.903
0.853
0.854
0.906
0.846
0.837
0.892
0.842
0.815
1688

Wheat Maize Other
Cereals
0.950 0.933 0.981
0.947 0.935 0.968
0.937 0.949 0.973
0.952 0.936 0.965
0.943 0.922 0.968
0.935 0.939 0.961
0.951 0.938 0.982
0.942 0.934 0.963
0.948 0.944 0.967
0.923 0.902 0.942
0.885 0.869 0.907
0.873 0.865 0.899
0.921 0.898 0.937
0.888 0.852 0.903
0.876 0.859 0.895
0.911 0.913 0.952
0.902 0.887 0.903
0.887 0.846 0.897
0.869 0.852 0.900
0.840 0.806 0.838
0.822 0.805 0.822
0.882 0.862 0.904
0.836 0.792 0.820
0.823 0.793 0.819
0.862 0.858 0.899
0.834 0.796 0.810
0.806 0.775 0.781
2380 1764 692

Fodder Perman. Other Woodl.
Crops Crops Cropl.
0.937 0.939 0.961 0.819
0.937 0.903 0.954 0.801
0.921 0.916 0.941 0.796
0.942 0.936 0.963 0.825
0.917 0.915 0.942 0.796
0.925 0.911 0.954 0.793
0.948 0.915 0.958 0.820
0.928 0.919 0.960 0.793
0.917 0.915 0.948 0.793
0.902 0.875 0.934 0.784
0.873 0.849 0.889 0.737
0.868 0.848 0.885 0.742
0.898 0.900 0.918 0.780
0.859 0.850 0.880 0.733
0.856 0.825 0.883 0.720
0.903 0.882 0.915 0.776
0.870 0.840 0.878 0.723
0.862 0.830 0.878 0.722
0.860 0.836 0.872 0.730
0.823 0.780 0.817 0.679
0.814 0.770 0.815 0.681
0.867 0.846 0.878 0.734
0.810 0.771 0.811 0.673
0.807 0.760 0.796 0.659
0.868 0.840 0.884 0.734
0.800 0.758 0.807 0.663
0.791 0.741 0.795 0.638
3120
1452 1368 5976

Grassl. Ot. Land Average
Cover
SBI
0.967 0.976
0.288
0.950 0.942
0.110
0.954 0.955
0.070
0.965 0.954
0.054
0.953 0.953
0.076
0.963 0.935
0.044
0.970 0.964
0.028
0.957 0.961
0.067
0.963 0.937
0.044
0.921 0.927
0.296
0.890 0.891
0.103
0.901 0.894
0.070
0.926 0.930
0.055
0.891 0.877
0.075
0.877 0.886
0.045
0.911 0.935
0.036
0.889 0.883
0.066
0.887 0.881
0.043
0.898 0.888
0.298
0.845 0.837
0.109
0.839 0.825
0.078
0.905 0.898
0.066
0.838 0.819
0.085
0.834 0.820
0.056
0.893 0.889
0.052
0.823 0.817
0.076
0.806 0.798
0.053
2368
1300

assurance of having respected the πi s, otherwise the HT estimator could provide extremely
biased results. From the CV (π̂i )s (9) reported in Table (3) it is clear that the standardization
of the distance matrix implied again as a result that the πs s are approximately constant even
if, with the exception of the LUCAS frame that lies on a regular grid, the error introduced
is quite higher than that observed in regular grid populations. The spatial distribution of the
population is thus an essential aspect to entail an efficient standardization of DU . In addition to the effects of γ and n already observed in regular grids it is possible to notice that a
combination of a small sample size, an high tuning parameter and the presence of clustering
in the coordinates of units in U may worryingly reduce the accuracy of the πi s. In such cases
it would be advisable to proceed to a partition of the frame population in zones with a more
homogeneous spatial distribution and then standardize the DU independently within each
stratum. It should also be considered that the CV (π̂i )s that may seem high actually did not
result in significant biases in the HT estimates proving its robustness to even medium-high
deviations in the πi s or in their estimated counterpart π̂i s. It is indeed necessary to remember that in this case we used only 10,000 sample replications instead of the 100,000 used
for regular grids, and probably they are not sufficient to achieve reliable estimates for the
inclusion probabilities.

Fast Selection of Spatially Balanced Samples

19

Looking at Tables (4), (5) and (6), it can be noticed that generally the HPWD gives encouraging results as it seems to handle any existing spatial data structure and effectively treat it to
locate units in the study region. Clustering of the population, presence of a spatial trend and
of spatial autocorrelation imply a clear direct effect on the reduction of estimation variance
even though their joint impact is obviously extremely moderate.
Especially a linear trend has been a valuable attribute to be exploited by our design even if
very high variance decreases are also found in the presence of spatial autocorrelation, while
the clustering of the population units slightly mitigate its capacity to spread the sample. Although it seems to be so sensitive to the occurrence of any of these properties it is also quite
robust in case of their absence, as it always has a RMSE much smaller to that obtained using
the SRS. Note that GRTS, SCPS and LPM have a similar behavior but with a lower gain in
efficiency in all these situations, in particular GRTS that shows to achieve samples that are
less spatially balanced than any other method.
The Meuse example confirms that HPWD is at least as efficient as PWD and that all the
trends found in artificial populations are also verified on this case study. The results of the
LUCAS data are also very encouraging. The gain in RMSE of the estimates compared with
SRS is remarkable reaching and exceeding, when n = 600 in some land cover classes, 25%.
In any experiment the PWD is slightly but systematically more efficient than SCPS and LPM
if used with γ ≥ 5 and sometimes a bit less efficient when used with lower exponents. The
HPWD proves to be an appreciable approximation of the PWD as it always has very similar
RMSEs and, in some rare cases, surprisingly better. The increase of n invariably entails an
improvement in the performance of PWD and HPWD with respect not only to SRS but also
to GRTS while the difference with SCPS and LPM seems to remain appreciably constant.
Table (7) reports the average CPU time in seconds of a 3,06 GHz Intel Core 2 Duo used by
each of the algorithms to select a sample for different populations and sample sizes. For each
method the time of both the developed version using R and C ++, when available (except for
the GRTS, both versions are always available), is reported. For SCPS and LPM, the R code
has been extracted from the available online material included in [26, 28]. The CPU time of
LPM are reported for both the suggested methods: LPM1 is slower but more accurate (used
to produce any results in this article) and LPM2 is faster but less accurate.
The extent to which the execution time matters, particularly if it highly depends on N and n,
is not secondary to the choice of the design. It is clear that among the examined procedures
the GRTS is the more computationally intensive method and that the HPWD is sensibly
quicker than the other distance based methods. The time spent to select a sample gradually
increases with n and only proportionally with N mainly because the number of iterations is
exactly n while the time used by SCPS and the two versions of the LPM increases with N 2 .
To be honest, however, this comparison is strongly influenced by the computation of DU
that in PWD and HPWD is excluded from the computing times, since, having this matrix to
be standardized in each application, its external calculation is required. In SCPS and LPM
the distances calculation is instead included in the sample selection algorithm.
Finally we can be confident that HPWD can be effectively applied without extensive difficulties even to spatial populations of large size and the only limit can be represented by the
amount of memory needed to store the distance matrix.

5 Conclusions
The use of space as a criterion for selecting sample units from a population is a solution that
promises significant developments in ecological, environmental, agricultural and forestry

20

F. Piersimoni, R. Benedetti

Table 7 Average computing time (in seconds) to select a sample based on 1,000 replications for each design
and for different population and sample sizes.

N
1000
1000
1000
1500
1500
1500
2000
2000
2000
2500
2500
2500
3000
3000
3000

R
C++
n GRTS SCPS LPM1 LPM2 PWD HPWD SCPS LPM1 LPM2 PWD HPWD
100 0.182 3.114 8.348 2.008 0.243 0.014 0.058 0.013 0.007 0.009 0.000
300 1.941 3.094 8.457 2.010 0.374 0.039 0.058 0.013 0.007 0.027 0.002
600 57.008 2.962 8.137 1.960 0.594 0.075 0.058 0.013 0.007 0.040 0.003
100 0.116 6.728 17.698 4.393 0.376 0.019 0.135 0.028 0.015 0.036 0.007
300 0.848 6.665 16.662 4.161 0.598 0.059 0.133 0.028 0.015 0.014 0.002
600 23.195 6.668 17.627 4.375 0.952 0.112 0.134 0.028 0.015 0.030 0.029
100 0.104 11.920 31.587 7.757 0.504 0.026 0.246 0.050 0.026 0.049 0.001
300 1.088 11.848 31.586 7.682 0.791 0.074 0.246 0.050 0.026 0.023 0.001
600 5.267 11.718 31.084 7.661 1.264 0.147 0.245 0.050 0.026 0.038 0.003
100 0.115 18.511 49.659 12.036 0.642 0.031 0.392 0.078 0.041 0.014 0.001
300 1.108 18.304 49.919 12.051 1.019 0.094 0.392 0.078 0.041 0.051 0.011
600 5.021 18.318 49.666 12.053 1.652 0.185 0.392 0.078 0.041 0.062 0.013
100 0.134 26.623 71.000 17.255 0.769 0.038 0.574 0.112 0.059 0.058 0.001
300 1.418 26.412 71.197 17.343 1.242 0.115 0.573 0.113 0.059 0.028 0.005
600 4.688 26.532 66.318 16.284 1.972 0.223 0.564 0.113 0.059 0.054 0.011

surveys. In addition to collecting data on phenomena that are impossible or particularly expensive to observe with an exhaustive direct observation of the frame population, this type
of surveys allows also to update existing lists frames (e.g. agricultural holdings), to use auxiliary information available only on a geographical basis (remotely sensed data), to facilitate
some aspects of quality controls and to better define some concepts and nomenclatures for
data dissemination (small area estimation).
Recent advances in geocoding that have led to a better understanding of the position of population units have thus extended the interest in these applications by the Institutes responsible
for the production of official statistics. Typical households and business surveys have been
based on archives that, in addition to many administrative data, include spatial information
on each single unit of the population. It also results that these institutes have to update the
survey techniques adopted so far. Among these, random sample selections have seen both
theoretical and practical developments in methods and algorithms, allowing us to conduct
surveys efficiently on spatial populations. Their motivation is still unclear and their framing
is still debated within two approaches, based on the model or on the design, formally very
different but that lead to logical and intuitive choices that can be reconciled. The selection
with probability proportional to the distance between sampling units seems to go in this direction with the use of the PWD method. One of the main limitations of this studies lies on
the fact that only the specification of homogeneous and isotropic spatial processes can be
considered by using the distance matrix as a summary of the population spatial features. The
population size, however, could be a big obstacle to be circumvented as this method is based
on a computationally intensive MCMC approximation of the distribution P (S ) representing
the sample design. Simultaneously increasing the amount of surveys on which to use these
methods and information details on the phenomena to be investigated, both in terms of spatial resolution and size, it has become essential to propose a much faster alternative to PWD,
the HPWD that can actually randomly select S in exactly n steps.
The results of the empirical studies show that the proposed HPWD design can compete in
precision with the best spatially balanced designs analyzed, while also allowing large reductions in the computational effort required. It represents an important alternative to such

Fast Selection of Spatially Balanced Samples

21

designs when the survey requires, as it is often the case in real situations, the use of a
very large spatial population. The results obtained seem to confirm that the simplicity of a
draw–by–draw scheme may also be sufficient to handle many spatial aspects of the frame
population and of the collected variables.
Some issues remain open for future research fundamentally linked to the theoretical derivation of πi s and πij s when the HPWD or PWD are used. Their knowledge would enable to
improve the study the theoretical properties of these designs and perhaps to derive more
flexible alternatives that can be better adapted to the treatment of spatial samples repeated
over time or to the selection of samples well–spread, stratified and simultaneously balanced
on some auxiliaries.

References
1. Arbia, G.: The use of GIS in spatial statistical surveys. International Statistical Review 61, 2, 339–359.
(1993)
2. Benedetti, R., Bee, M., Espa, G., Piersimoni, F. (eds): Agricultural Survey Methods. Chicester: Wiley.
(2010)
3. Benedetti, R., Espa, G., Taufer, E.: Model-based variance estimation in non-measurable spatial designs,
Journal of Statistical Planning and Inference 181, 52–61. (2017)
4. Benedetti, R., Palma, D.: Optimal sampling designs for dependent spatial units. Environmetrics 6, 101–
114. (1995)
5. Benedetti, R., Piersimoni, F.: A spatially balanced design with probability function proportional to the
within sample distance. Biometrical Journal doi:10.1111/insr.12216. (2017)
6. Benedetti, R., Piersimoni, F., Postiglione, P.: Sampling spatial units for agricultural surveys. Advances in
Spatial Science Series. Berlin, Heidelberg: Springer-Verlag. (2015)
7. Benedetti, R., Piersimoni, F., Postiglione, P.: Spatially balanced sampling: a review and a reappraisal.
International Statistical Review doi:10.1111/insr.12216. (2017)
8. Breidt, F.J., Chauvet, G.: Penalized balanced sampling. Biometrika 99, 945–958. (2012)
9. Chambers, R.L., Clark, R.G.: An introduction to model–based survey sampling with applications. Oxford:
Oxford University Press. (2012)
10. Chauvet, G., Tillè, Y.: A fast algorithm of balanced sampling. Computational Statistics 21, 53–62. (2006)
11. Cressie, N.A.C.: Statistics for spatial data. 2nd edition. New York: Wiley. (1993)
12. Cressie, N.A.C., Chambers, R. L..: Comment on Article by Ferreira and Gamerman. Bayesian Analysis
10, 3, 741–748. (2015)
13. Cressie, N.A.C., Johannesson G.: Fixed rank kriging for very large spatial data sets. Journal of the Royal
Statistical Society, Series B 70, 209–226. (2008)
14. Cordy, C.B.: An extension of the HorvitzThompson theorem to point sampling from a continuous universe. Statististics and Probability Letters 18, 353–362. (1993)
15. Cromley, E., McLafferty, S.: GIS and Public Health. New York: Guilford Press. (2002)
16. Deville, J.C., Tillè, Y.: Unequal probability sampling without replacement through a splitting method.
Biometrika 85, 89–101. (1998)
17. Deville, J.C., Tillè, Y.: Efficient balanced sampling: the cube method. Biometrika 91, 893–912. (2004)
18. Diggle, P., Menezes, R., Su, T.: Geostatistical inference under preferential sampling. Journal of the Royal
Statistical Society, Series C (Applied Statistics) 59, 191–232. (2010)
19. Diggle, P.J., Ribeiro, P.J.: Model-based geostatistics. New York: Springer-Verlag. (2007)
20. Fattorini, L.: Applying the Horvitz–Thompson criterion in complex designs: A computer-intensive perspective for estimating inclusion probabilities. Biometrika 93, 269-278. (2006)
21. Fattorini, L.: An adaptive algorithm for estimating inclusion probabilities and performing the Horvitz–
Thompson criterion in complex designs. Computational Statistic 24, 623-639. (2009)
22. Ferreira, G., Gamerman, D.: Optimal design in geostatistics under preferential sampling. Bayesian Analysis 10, 3, 711–735. (2015)
23. Gaetan, C., Guyon, X.: Spatial statistics and modeling. Springer series in statistics. New York: SpringerVerlag. (2010)
24. Gallego, J., Delincé, J.: The European land use and cover area-frame statistical survey, In: Agricultural
survey methods. Benedetti, R., Bee, M., Espa, G. and Piersimoni, F. (eds), John Wiley and Sons Ltd, Chichester, UK, 151–168. (2010)

22

F. Piersimoni, R. Benedetti

25. Gelfand, A.E., Diggle, P.J., Fuentes, M., Guttorp, P.: Handbook of spatial statistics. Handbooks of modern statistical methods. Boca Raton: Chapman and Hall/CRC. (2010)
26. Grafström, A.: Spatially correlated Poisson sampling. Journal of Statistical Planning and Inference 142,
139–147. (2012)
27. Grafström, A., Lisic, J.: BalancedSampling: balanced and spatially balanced sampling. R package
version 1.5.2. https://CRAN.R-project.org/package=BalancedSampling. Accessed July 10, 2017. (2017)
28. Grafström, A., Lundström, N.L.P., Schelin, L.: Spatially balanced sampling through the pivotal method.
Biometrics 68, 514–520. (2012)
29. Grafström, A., Tillè, Y.: Doubly balanced spatial sampling with spreading and restitution of auxiliary
totals. Environmetrics 24, 120–131. (2013)
30. Gregoire, T.G., Valentine, H.T.: Sampling strategies for natural resources and the environment. Applied
environmental statistics. London: Chapman and Hall/CRC (2008)
31. Hainy, M., Müller, W.G., Wagner, H.: Likelihood-free simulation-based optimal design with an application to spatial extremes. Stochastic Environmental Research and Risk Assessment 30, 2, 481–492. (2016)
32. Hainy, M., Müller, W.G., Wynn, H. P.: Learning functions and approximate bayesian computation design:
ABCD. Entropy 16, 8, 4353–4374. (2014)
33. Hedayat, A., Rao, C.R., Stufken, J.: Sampling designs excluding contiguous units. Journal of Statistical
Planning and Inference 19, 159–170. (1988)
34. Isaki, C.T., Fuller, W.A.: Survey design under the regression superpopulation model. Journal of the
American Statistical Association 77, 377, 89–96. (1982)
35. Johannesson, G., Cressie, N.A.C.: Finding large-scale spatial trends in massive, global, environmental
datasets. Environmetrics 15, 1, 1–44. (2004)
36. Johnson, C.R., Masson, R.D., Trosset, M.W.: On the diagonal scaling of Euclidean distance matrices to
doubly stochastic matrices. Linear Algebra and its Applications 397, 253-264. (2005)
37. Jones, C.B.: Geographical Information Systems and Computer Cartography. Harlow, Essex, UK: Addison Wesley Longman. (1997)
38. Kincaid, T.M., Olsen, A.R.: spsurvey: spatial survey design and analysis. R package version 3.3. (2016)
39. Knottnerus, P.: Sample survey theory: some Pythagorean perspectives. New York: Springer-Verlag.
(2003)
40. LeSage, J.P., Pace, R.K.: Introduction to spatial econometrics. Statistics, textbooks and monographs.
Boca Raton: Chapman and Hall/CRC. (2009)
41. Mandallaz, D.: Sampling techniques for forest inventories. Applied environmental statistics. Boca Raton:
Chapman and Hall/CRC. (2008).
42. Möller J.: Spatial Statistics and Computational Method. New York: Springer-Verlag. (2003)
43. Müller, P.: Simulation-Based optimal design. Bayesian Statistics 6, 459–474. (1999)
44. Müller, W.G.: Collecting spatial data: optimum design of experiments for random fields. Berlin:
Springer-Verlag, 3rd edition (2007)
45. Müller, W.G., Pronzato, L., Rendas, J., Waldl, H.: Efficient prediction designs for random fields. Applied
Stochastic Models in Business and Industry 31, 2, 178–194. (2015)
46. Müller, W.G., Stehlı́k, M.: Compound optimal spatial designs, Environmetrics 21, 3-4, 354–364. (2010)
47. R Core Team: R: a language and environment for statistical computing. R Foundation for Statistical
Computing. Vienna, Austria. Available from:https://www.R-project.org/. Accessed March 10, 2017. (2016)
48. Rogerson, P., Delmelle, E.: Optimal sampling design for variables with varying spatial importance. Geographical Analysis 36, 177–194. (2004)
49. Särndal, C.E., Swensson, B. and Wretman, J.: Model Assisted Survey Sampling. New York: SpringerVerlag. (1992)
50. Schabenberger, O., Gotway, C.A.: Statistical Methods for Spatial Data Analysis. Texts in Statistical
Science Series, Boca Raton: Chapman and Hall/CRC. (2004)
51. Scheiner, S.M., Gurevitch, J.: Design and Analysis of Ecological Experiments. 2nd edition. Oxford, UK:
Oxford University Press. (2001)
52. Stevens, Jr.D.L., Olsen, A.R.: Variance estimation for spatially balanced samples of environmental resources. Environmetrics 14, 593–610. (2003)
53. Stevens, Jr.D.L., Olsen, A.R.: Spatially balanced sampling of natural resources. Journal of the American
Statistical Association 99, 262–278. (2004)
54. Tillé, Y.: Sampling Algorithms. New York: Springer-Verlag. (2006)
55. Tillé, Y.: Ten years of balanced sampling with the CUBE method: An appraisal. Survey Methodology 2,
215–226. (2011)
56. Tillé, Y., Wilhelm, M.: Probability Sampling Designs: Principles for Choice of Design and Balancing.
Statistical Science 32, 2, 176–189. (2017)
57. Tzeng, S., Huang, H.C., Cressie, N.A.C.: A Fast, Optimal Spatial–Prediction Method for Massive
Datasets. Journal of the American Statistical Association 100, 472, 1343–1357. (2005)

Fast Selection of Spatially Balanced Samples

23

58. Valliant, R., Dorfman. A.H., Royall, R.M.: Finite population sampling and inference: a prediction approach. New York: John Wiley and Sons. (2000)
59. Waagepetersen, R.: An estimating function approach to inference for inhomogeneous Neyman-Scott
processes. Biometrics 63, 252–258. (2007)
60. Wang, J.F., Jiang,C.S., Hu, M.G., Cao, Z.D., Guo, Y.S., Li, L.F., Liu, T.J., Meng, B.: Design-based spatial
sampling: Theory and implementation. Environmental Modelling and Software 40, 280–288. (2013)
61. Wang, J.F., Stein, A., Gao, B.B., Ge, Y.: A review of spatial sampling. Spatial Statistics 2, 1–14. (2012)
62. Webster, R., Oliver, M.A.: Geostatistics for Environmental Scientists. 2nd Edition, Statistics in Practice.
Chichester: John Wiley and Sons. (2007)

