Getting Started with sondage

Without replacement

The simplest case: draw n units from a population of size N with equal probabilities.

library(sondage)
s <- equal_prob_wor(N = 50, n = 5)
s
#> Equal prob WOR [srs] (n=5, N=50): 4 39 1 34 23

The returned object stores the sampled unit indices and the inclusion probabilities:

s$sample
#> [1]  4 39  1 34 23
s$pik[1:10]  # all equal to n/N = 0.1
#>  [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

Alternative methods are available. Systematic sampling provides implicit stratification based on the ordering of units:

s_sys <- equal_prob_wor(N = 50, n = 5, method = "systematic")
s_sys
#> Equal prob WOR [systematic] (n=5, N=50): 7 17 27 37 47

Bernoulli sampling selects each unit independently with probability p = n/N, producing a random sample size. The parameter n is the expected sample size not the realized sample size.

s_ber <- equal_prob_wor(N = 50, n = 5, method = "bernoulli")
s_ber
#> Equal prob WOR [bernoulli] (expected n=5, realized n=2, N=50): 10 13
length(s_ber$sample)  # may differ from 5
#> [1] 2

With replacement

s_wr <- equal_prob_wr(N = 50, n = 5)
s_wr
#> Equal prob WR [srs] (n=5, N=50): 39 42 6 24 32

With-replacement designs track hits (how many times each unit was selected) instead of inclusion probabilities:

s_wr$hits  # length-N vector of selection counts
#>  [1] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0
#> [39] 1 0 0 1 0 0 0 0 0 0 0 0

Computing design quantities

Start from a vector of positive size measures and compute inclusion probabilities (WOR) or expected hits (WR):

x <- c(10, 20, 5, 40, 25)
n <- 3

# Without replacement: inclusion probabilities (capped at 1)
pik <- inclusion_prob(x, n)
pik
#> [1] 0.3333333 0.6666667 0.1666667 1.0000000 0.8333333
sum(pik)  # equals n
#> [1] 3

# With replacement: expected hits (can exceed 1)
eh <- expected_hits(x, n)
eh
#> [1] 0.30 0.60 0.15 1.20 0.75
sum(eh)   # equals n
#> [1] 3

Without replacement

Pass the inclusion probability vector to unequal_prob_wor(). The default method is Conditional Poisson Sampling (CPS, Chen et al. 1994), which has maximum entropy among fixed-size designs:

s <- unequal_prob_wor(pik, method = "cps")
s
#> Unequal prob WOR [cps] (n=3, N=5): 1 4 5
s$sample
#> [1] 1 4 5
inclusion_prob(s)  # extract the stored design-defining target pik
#> [1] 0.3333333 0.6666667 0.1666667 1.0000000 0.8333333

Other fixed-size methods include "brewer" (draw-by-draw, order-invariant), "sampford" (Sampford 1967), "systematic" (O(N), implicit stratification), "sps" (Sequential Poisson sampling, Ohlsson 1998, supports PRN), and "pareto" (Pareto sampling, Rosén 1997, supports PRN). The "poisson" method produces a random sample size.

For sps and pareto, inclusion_prob(s) still returns the stored target pik vector used to define the design. Those methods do not, in general, provide exact finite-population first-order inclusion probabilities equal to that target.

With replacement

Pass the expected hits vector to unequal_prob_wr(). Chromy (1979)’s method is the default:

s_wr <- unequal_prob_wr(eh, method = "chromy")
s_wr
#> Unequal prob WR [chromy] (n=3, N=5): 2 4 5
s_wr$hits  # realized selection counts
#> [1] 0 1 0 1 1

The "multinomial" method draws each selection independently.

The "systematic_mr" method is systematic PPS with minimum replacement. It walks the cumulated expected hits from one random start, in the order of the frame. A unit larger than the sampling interval is hit once per interval it spans, where unequal_prob_wor(method = "systematic") takes it once, with certainty. Here the last unit expects 1.875 hits, so it receives one or two:

eh_big <- expected_hits(c(10, 20, 5, 40, 125), n)
eh_big
#> [1] 0.150 0.300 0.075 0.600 1.875
unequal_prob_wr(eh_big, method = "systematic_mr")$hits
#> [1] 0 0 1 0 2

Stratified balanced sampling

When strata are known, the stratified cube method (Chauvet 2009) preserves within-stratum sample sizes while balancing on auxiliary variables. Exact preservation requires sum(pik) within each stratum to be close to an integer. When this condition is not met, within-stratum and total sample sizes will vary around their targets and the design is treated as random-size (fixed_size = FALSE):

N <- 20
pik_strat <- rep(0.4, N)
x_strat <- matrix(as.double(1:N), ncol = 1)
strata <- rep(1:4, each = 5)

s_strat <- balanced_wor(pik_strat, aux = x_strat, strata = strata)
s_strat
#> Balanced WOR [cube] (n=8, N=20): 1 5 7 10 12 13 17 20

# Check per-stratum sample sizes
table(strata[s_strat$sample])
#> 
#> 1 2 3 4 
#> 2 2 2 2

Joint inclusion probabilities for the cube use the high-entropy approximation. It depends only on pik, so it ignores strata and balancing constraints, and it can give a positive probability to a pair that a stratified design excludes.

Balanced sampling with inequalities

The cube method also accepts linear inequality constraints through bounds = list(B, lower, upper) (Tripet and Tillé 2026). These constrain raw realized-sample sums, \(\mathtt{lower}_j \leq \sum_{k \in S} B_{kj} \leq \mathtt{upper}_j\), while preserving the supplied first-order inclusion probabilities. A common use is to keep category counts between the integers immediately below and above their expectations:

pik_control <- rep(0.5, 12)
group <- rep(c("a", "b", "c"), each = 4)
B <- sapply(unique(group), function(g) as.double(group == g))
expected_counts <- colSums(B * pik_control)

s_control <- balanced_wor(
  pik_control,
  bounds = list(
    B = B,
    lower = floor(expected_counts),
    upper = ceiling(expected_counts)
  )
)
colSums(B[s_control$sample, , drop = FALSE])
#> a b c 
#> 2 2 2

Spatially balanced sampling

It’s possible to select well-spread, spatially balanced samples using lpm2 (local pivotal method 2, Grafström et al. 2012) and scps (spatially correlated Poisson sampling, Grafström 2012).

coordinates <- as.matrix(expand.grid(x = 1:4, y = 1:4))
pik_spatial <- rep(0.25, nrow(coordinates))

s_lpm2 <- balanced_wor(
  pik_spatial,
  spread = coordinates,
  method = "lpm2"
)
s_scps <- balanced_wor(
  pik_spatial,
  spread = coordinates,
  method = "scps"
)
s_lpm2$sample
#> [1]  1  6 12 15
s_scps$sample
#> [1]  2  4 11 14

These built-in methods are spread-only: they do not accept aux, strata, or bounds. They preserve the supplied first-order inclusion probabilities, but intentionally produce low-entropy designs, so joint_inclusion_prob() and sampling_cov() are not available for them.

With permanent random numbers, scps visits the units in row order and selects unit \(k\) when \(u_k\) falls below its current conditional probability (Bondesson and Thorburn 2008; Grafström and Matei 2018). Two surveys drawn with prn = u and prn = 1 - u are then negatively coordinated. Matei et al. (2023) spread on a measure of response burden instead of coordinates, so that highly burdened units rarely enter the same sample, and rarely enter both.

set.seed(42)
burdened <- as.double(seq_len(16) %in% c(2, 7, 11))
u <- runif(16)
s1 <- balanced_wor(pik_spatial, spread = burdened, method = "scps", prn = u)
s2 <- balanced_wor(pik_spatial, spread = burdened, method = "scps", prn = 1 - u)
intersect(s1$sample, s2$sample)
#> integer(0)

Joint inclusion probabilities (WOR)

pik <- inclusion_prob(c(10, 20, 5, 40, 25), n = 3)
s <- unequal_prob_wor(pik, method = "cps")

pikl <- joint_inclusion_prob(s)
pikl
#>            [,1]       [,2]       [,3]      [,4]      [,5]
#> [1,] 0.33333333 0.10026673 0.01935742 0.3333333 0.2137092
#> [2,] 0.10026673 0.66666667 0.04704252 0.6666667 0.5193574
#> [3,] 0.01935742 0.04704252 0.16666667 0.1666667 0.1002667
#> [4,] 0.33333333 0.66666667 0.16666667 1.0000000 0.8333333
#> [5,] 0.21370918 0.51935742 0.10026673 0.8333333 0.8333333

The diagonal holds the first-order probabilities \(\pi_i\) and the off-diagonal entries are \(\pi_{ij} = P(i \in S \text{ and } j \in S)\). For sps and pareto, the approximation is built from the stored target pik vector, not from unknown exact finite-population first-order inclusion probabilities.

Joint expected hits (WR)

eh <- expected_hits(c(10, 20, 5, 40, 25), n = 3)
s_wr <- unequal_prob_wr(eh, method = "chromy")

joint_expected_hits(s_wr)
#>        [,1]   [,2]   [,3]   [,4]   [,5]
#> [1,] 0.2968 0.1126 0.0332 0.3504 0.0974
#> [2,] 0.1126 0.6011 0.0103 0.6382 0.4411
#> [3,] 0.0332 0.0103 0.1549 0.1630 0.1033
#> [4,] 0.3504 0.6382 0.1630 1.6063 0.8484
#> [5,] 0.0974 0.4411 0.1033 0.8484 0.7451

Submatrices for large populations

By default, joint_inclusion_prob() and joint_expected_hits() return the full N x N matrix, which is prohibitive when N > 10,000 (~0.8 GB). In practice, variance estimators only need the submatrix for sampled units. The sampled_only = TRUE argument returns this n x n submatrix directly:

pik_large <- inclusion_prob(runif(500), n = 10)
s_large <- unequal_prob_wor(pik_large, method = "brewer")

# Full 500 x 500 matrix
pikl_full <- joint_inclusion_prob(s_large)
dim(pikl_full)
#> [1] 500 500

# Only the 10 x 10 submatrix for sampled units
pikl_sub <- joint_inclusion_prob(s_large, sampled_only = TRUE)
dim(pikl_sub)
#> [1] 10 10

All methods compute the submatrix directly without allocating the N x N matrix, so sampled-only queries reach well beyond the population size that a full matrix allows. Memory and computing time still depend on the method, N and n, and the sampled output is limited to 10,000 units. The sampling_cov() generic also accepts sampled_only:

delta_sub <- sampling_cov(s_large, sampled_only = TRUE)
dim(delta_sub)
#> [1] 10 10

Sampling covariance

The sampling_cov() generic computes \(\Delta_{ij} = \pi_{ij} - \pi_i \pi_j\) for WOR designs (or the analogous expression for WR designs):

sampling_cov(s)
#>             [,1]        [,2]        [,3] [,4]        [,5]
#> [1,]  0.22222222 -0.12195549 -0.03619814    0 -0.06406859
#> [2,] -0.12195549  0.22222222 -0.06406859    0 -0.03619814
#> [3,] -0.03619814 -0.06406859  0.13888889    0 -0.03862216
#> [4,]  0.00000000  0.00000000  0.00000000    0  0.00000000
#> [5,] -0.06406859 -0.03619814 -0.03862216    0  0.13888889

With weighted = TRUE, it returns the Sen-Yates-Grundy check quantities \(1 - \pi_i \pi_j / \pi_{ij}\) (Sen 1953; Yates and Grundy 1953). For well-behaved WOR designs, these should be non-positive off-diagonal:

sampling_cov(s, weighted = TRUE)
#>            [,1]        [,2]       [,3] [,4]        [,5]
#> [1,]  0.6666667 -1.21631058 -1.8699881    0 -0.29979336
#> [2,] -1.2163106  0.33333333 -1.3619296    0 -0.06969793
#> [3,] -1.8699881 -1.36192957  0.8333333    0 -0.38519411
#> [4,]  0.0000000  0.00000000  0.0000000    0  0.00000000
#> [5,] -0.2997934 -0.06969793 -0.3851941    0  0.16666667

Notes on specific methods

Choosing a method

• CPS: Use when exactness matters more than speed. It is the maximum-entropy fixed-size unequal-probability design, with exact first- and second-order inclusion probabilities.
• Brewer: Use when you want exact first-order inclusion probabilities with lower computational cost than cps, and can work with approximate second-order quantities in sondage.
• Sampford: Use for fixed-size PPS with exact first- and second-order inclusion probabilities and fast typical draws.
• Systematic: Use when very fast sampling and ordering or implicit stratification are central, and structural zeros in some joint inclusion probabilities are acceptable (low-entropy).
• Systematic PPS with minimum replacement: Use when a large unit must be able to receive several hits, for instance to keep a fixed second-stage take per hit.
• Poisson: Use when a random sample size is acceptable and independent selection is desirable.
• SPS / Pareto: Use as fast high-entropy order-sampling alternatives when approximate first- and second-order quantities are acceptable.
• Cube: Use when balancing on auxiliary variables is important.
• LPM2 / SCPS: Use when spatial spread is the priority and a local-neighborhood variance strategy is available.

Choosing a variance strategy

• SYG with exact \(\pi_{ij}\): CPS, Sampford and SRS WOR provide exact joint probabilities. For these fixed-size designs, the Sen-Yates-Grundy variance estimator (Sen 1953; Yates and Grundy 1953) is unbiased when every joint inclusion probability is positive, which excludes n = 1.
• SYG with approximate \(\pi_{ij}\): For Brewer, SPS, Pareto, and cube, joint probabilities use the high-entropy approximation (Brewer and Donadio 2003, eq. 18). This guarantees symmetry, non-negativity, and \(\pi_{ij} \le \min(\pi_i, \pi_j)\), but does not satisfy the fixed-size marginal identity \(\sum_{j \neq i}\pi_{ij} = (n-1)\pi_i\) exactly. The marginal defect is typically small but can be non-trivial for highly skewed \(\pi_k\) vectors when \(n\) is small (e.g., \(n = 2\)). A warning is issued when it exceeds 5% of \(n\). For sps and pareto, this approximation is built from the stored target pik vector rather than the unknown exact finite-population first-order probabilities. For the cube, it also ignores strata and balancing constraints, as noted above. Use method = "cps" when exact second-order structure is required.
• Independent designs (Poisson, Bernoulli): The covariance \(\Delta_{ij} = 0\) for \(i \neq j\), so the variance simplifies to \(\text{Var}(\hat{Y}_{HT}) = \sum_i (1 - \pi_i)\, y_i^2 / \pi_i\).
• Systematic designs with \(\pi_{ij} = 0\): The SYG estimator is undefined. Use successive-differences, Hartley-Rao approximations, or replication methods.
• With-replacement designs (Chromy, multinomial, systematic PPS with minimum replacement): Use the Hansen and Hurwitz (1943) estimator, which is design-unbiased for multinomial draws and an approximation, often conservative, for the two minimum-replacement methods, or the generalized Horvitz-Thompson framework of Chromy (2009) which unifies WOR and WR variance estimation via sampling_cov().