Balanced Sampling Without Replacement

Draws a balanced sample using the cube method (Deville & Tillé, 2004), or a spatially balanced (well-spread) sample using the local pivotal method (Grafström, Lundström & Schelin, 2012) or spatially correlated Poisson sampling (Grafström, 2012). A balanced sample satisfies (approximately) the balancing equations $\sum_{k \in S} x_k / \pi_k \approx \sum_{k \in U} x_k$ for each auxiliary variable $x$; a well-spread sample selects units that are far apart in the space spanned by the spreading variables.

Usage

balanced_wor(
  pik,
  aux = NULL,
  strata = NULL,
  spread = NULL,
  bounds = NULL,
  method = c("cube", "lpm2", "scps"),
  nrep = 1L,
  ...
)

Arguments

pik

A numeric vector of inclusion probabilities (length N). sum(pik) must be an integer to floating-point accuracy; see unequal_prob_wor() for the exact-0/1 handling of boundary values. The target sample size, sum(pik), must be at least 1.

aux

An optional numeric matrix (N x p) of auxiliary balancing variables. Each column defines a balancing constraint. The sample size constraint is always included automatically; aux specifies additional variables to balance on. When NULL, only the sample size is balanced (equivalent to an unbalanced fixed-size design).

strata

An optional integer vector (length N) of stratum indicators (positive integers). Uses the stratified cube method (Chauvet & Tillé, 2006; Chauvet, 2009) to preserve within-stratum sample sizes while balancing on aux. Requires sum(pik) within each stratum to be close to an integer for exact sizes. If not, a warning is issued and fixed_size is set to FALSE.

spread

An optional numeric matrix (N x d) of spatial coordinates (or other spreading variables) for well-spread, spatially balanced sampling. Required by the built-in "lpm2" and "scps" methods, and supported by methods registered via register_method() with supports_spread = TRUE; the built-in "cube" method does not use it and will error. A non-matrix vector is treated as a single spreading variable.

bounds

An optional list with elements B, lower, and upper describing linear inequality constraints on the realized sample (Tripet & Tillé, 2026): B is a numeric matrix (N x q) whose columns are constraint variables, and the sample s is drawn so that $lower_j \le \sum_{k \in s} B_{kj} \le upper_j$ for every constraint j. -Inf / Inf entries make a constraint one-sided; lower[j] == upper[j] enforces an exact equality. The starting probabilities must be feasible: lower <= colSums(B * pik) <= upper. Only supported by the built-in "cube" method. See Inequality constraints below.

method

The sampling method. "cube" (the default) balances on aux; "lpm2" (local pivotal method 2) and "scps" (spatially correlated Poisson sampling) spread on spread; or the name of a balanced method added via register_method().

nrep

Number of replicate samples (default 1). When nrep > 1, $sample holds a matrix (n x nrep) for fixed-size designs, or a list of integer vectors when within-stratum sizes are not exact.

...

Additional arguments passed to methods:

eps: Boundary tolerance (default 1e-10): decides when an updated working probability has numerically reached 0 or 1 (during the cube flight phase, or during the pivotal steps of "lpm2" or "scps"). It never reclassifies the input; supplied pik inside (0, eps] or [1 - eps, 1) are rejected. Only pik of exactly 0 or exactly 1 enter the algorithm as already-resolved units.
condition_aux: Logical; if TRUE, pre-conditions aux by weighted centering/scaling and QR-pivot rank pruning to improve numerical stability with ill-conditioned or collinear auxiliary variables (default FALSE).
qr_tol: Tolerance for QR rank detection when condition_aux = TRUE; must be one finite, non-negative number (default sqrt(.Machine$double.eps)).

Value

An object of class c("balanced", "unequal_prob", "wor", "sondage_sample"). When nrep = 1, $sample is an integer vector of selected unit indices. When nrep > 1, $sample is a matrix (n x nrep) for fixed-size designs, or a list of integer vectors when fixed_size is FALSE (e.g., stratified with non-integer per-stratum sizes).

Details

The cube method proceeds in two phases:

Flight phase: Probabilities are moved toward 0 or 1 while maintaining all balancing constraints. Each step resolves at least one unit. Terminates when fewer than p+1 undecided units remain.
Landing phase: Remaining undecided units are resolved by progressively relaxing balancing constraints, starting from the last column of aux. Users should order auxiliary variables by importance (most important first).

The sample size constraint is always placed first (never relaxed during landing). For stratified designs, within-stratum size constraints are also placed first.

Joint inclusion probabilities are approximated via the high-entropy approximation (Brewer & Donadio, 2003), which is appropriate since the cube produces a near-maximum-entropy design.

Inequality constraints

bounds implements the cube method with inequality constraints of Tripet & Tillé (2026). During the flight phase, steps are capped so every constraint stays feasible, and a constraint whose slack reaches zero becomes an equality from then on. The inclusion probabilities are respected exactly (E(s) = pik), unlike rejective procedures.

B applies to the realized sample directly (counts / raw sums); it is not divided by pik. To bound a Horvitz-Thompson estimator, pass x / pik as the constraint column.

The main application is controlled selection (Goodman & Kish, 1950): bounding category counts to the integers adjacent to their expectation. With indicator columns B and S <- colSums(B * pik), use lower = floor(S), upper = ceiling(S). Categories may overlap (e.g. row and column margins of a two-way control table, as in NAEP-style designs).

Integer-valued bound systems on partitions or two-way margins are satisfied exactly. For structures with no exact integer solution, some three-way controlled rounding problems, or continuous-valued constraints that end the flight phase tight against a boundary, bounds that provably block the landing phase are relaxed one at a time, with a warning; E(s) = pik still holds. In other words, the bounds are guaranteed whenever no relaxation warning is raised.

The high-entropy approximation used by joint_inclusion_prob() is less accurate under tight bounds, which distort the design away from maximum entropy; Tripet & Tillé (2026) recommend Monte Carlo estimation of joint inclusion probabilities in that case.

Spatially balanced sampling (`lpm2` and `scps`)

method = "lpm2" implements the local pivotal method 2 of Grafström, Lundström & Schelin (2012). Repeatedly, a randomly chosen undecided unit and its nearest undecided neighbour in the spread space compete in a pivotal step (Deville & Tillé, 1998) that resolves at least one of them to 0 or 1 while preserving the inclusion probabilities exactly (E(s) = pik). Nearby units thereby tend to exclude each other, spreading the sample over the population. Nearest-neighbour distance ties (common on gridded coordinates) are broken uniformly at random.

Spreading variables should be on comparable scales, since nearness is plain Euclidean distance in the spread columns; rescale them (e.g. with scale()) when they are not. "lpm2" and "scps" use spread only: they do not accept aux, strata, or bounds. To exactly balance covariate totals and spread spatially, register a method that supports both aux and spread (for example, a local-cube implementation). Running "cube" with coordinates in aux balances their totals but does not itself enforce spatial spread.

method = "scps" implements spatially correlated Poisson sampling with the maximal-weight strategy of Grafström (2012). At each step, a randomly chosen undecided unit is accepted or rejected using its current conditional probability. Its probability displacement is then distributed to the nearest undecided units, subject to feasibility bounds that keep every working probability in $[0, 1]$. Equal-distance units share weight as evenly as their bounds allow. Random selection of the step unit avoids dependence on input row order.

Both spatial methods deliberately drive joint inclusion probabilities of nearby units toward zero, so the design is not high entropy and no joint-probability approximation is provided: joint_inclusion_prob() errors for these designs. Variance for well-spread samples is usually estimated with local-neighbourhood estimators (e.g. Grafström & Schelin, 2014).

Both methods are spread-only: neither exactly balances the totals of aux. They are dispatched by balanced_wor() because “spatially balanced sampling” is the standard name for well-spread fixed-size designs, and because the same interface accommodates local-cube methods that combine exact balancing with spread. Capability metadata keeps the distinction explicit: both report supports_aux = FALSE and supports_spread = TRUE through method_spec().

LPM2 costs O(N^2 * d) time per draw. SCPS uses weighted quickselect to find the distance at which its maximal weights sum to one, avoiding a full sort of the remaining units at each step. Its expected cost is also O(N^2 * d). Both implementations use O(N) workspace and store no distance matrix; SCPS sorts only equal-distance cutoff groups to share their weight fairly.

References

Deville, J.C. and Tillé, Y. (1998). Unequal probability sampling without replacement through a splitting method. Biometrika, 85(1), 89-101.

Deville, J.C. and Tillé, Y. (2004). Efficient balanced sampling: the cube method. Biometrika, 91(4), 893-912.

Chauvet, G. and Tillé, Y. (2006). A fast algorithm for balanced sampling. Computational Statistics, 21(1), 53-62.

Chauvet, G. (2009). Stratified balanced sampling. Survey Methodology, 35, 115-119.

Grafström, A., Lundström, N.L.P. and Schelin, L. (2012). Spatially balanced sampling through the pivotal method. Biometrics, 68(2), 514-520. doi:10.1111/j.1541-0420.2011.01699.x

Grafström, A. (2012). Spatially correlated Poisson sampling. Journal of Statistical Planning and Inference, 142(1), 139-147. doi:10.1016/j.jspi.2011.07.003

Grafström, A. and Schelin, L. (2014). How to select representative samples. Scandinavian Journal of Statistics, 41(2), 277-290. doi:10.1111/sjos.12016

Tripet, A. and Tillé, Y. (2026). Balanced sampling with inequalities: application to category bounding, matrix rounding, and spread sampling. Journal of the American Statistical Association, 121(553), 796-806. doi:10.1080/01621459.2025.2550667

Examples

# Unequal probability balanced sample
pik <- c(0.3, 0.6, 0.2, 0.4, 0.5)
x <- matrix(c(10, 20, 15, 25, 30))
set.seed(1)
s <- balanced_wor(pik, aux = x)
s$sample
#> [1] 2 5

# Check balancing: HT estimate of aux totals vs population totals
colSums(x[s$sample, , drop = FALSE] / pik[s$sample]) - colSums(x)
#> [1] -6.666667

# Stratified balanced sample
N <- 20
pik <- rep(0.4, N)
x <- matrix(as.double(1:N), ncol = 1)
strata <- rep(1:4, each = 5)
set.seed(1)
s <- balanced_wor(pik, aux = x, strata = strata)
s$sample
#> [1]  1  2  6  8 13 15 16 19

# Controlled selection: bound category counts to the integers
# adjacent to their expectation (floor/ceil)
pik <- rep(0.5, 12)  # n = 6
groups <- rep(c("a", "b", "c"), each = 4)
B <- sapply(unique(groups), function(g) as.double(groups == g))
S <- colSums(B * pik)  # 2, 2, 2 per category
set.seed(1)
s <- balanced_wor(
  pik,
  bounds = list(B = B, lower = floor(S), upper = ceiling(S))
)
table(groups[s$sample])  # exactly 2 per category
#> 
#> a b c 
#> 2 2 2 

# Spatially balanced (well-spread) sample: local pivotal method 2
N <- 100
coords <- cbind(runif(N), runif(N))
pik <- rep(0.1, N)
set.seed(1)
s <- balanced_wor(pik, spread = coords, method = "lpm2")
s$sample
#>  [1]  2  5 15 21 38 40 46 49 81 85

# Spatially correlated Poisson sampling uses the same spread contract
set.seed(1)
s_scps <- balanced_wor(pik, spread = coords, method = "scps")
s_scps$sample
#>  [1]  1 10 12 17 23 27 61 75 90 97