--- title: "Getting Started with sondage" output: rmarkdown::html_vignette bibliography: references.bib vignette: > %\VignetteIndexEntry{Getting Started with sondage} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup, echo = FALSE} library(sondage) set.seed(1) ``` ## Overview `sondage` provides fast implementations of survey sampling algorithms for finite populations. Every sampling function returns a *design object* that carries the selected sample alongside the design metadata (inclusion probabilities or expected hits, sample size, method). A set of generics then operates on these objects to compute quantities needed for variance estimation. For the full object contract and field definitions, see `?sondage_sample`. For without-replacement designs, the stored `pik` vector is the design-defining target inclusion probability vector. For methods with exact first-order guarantees, this equals the true first-order inclusion probabilities. For order-sampling methods such as `sps` and `pareto`, the stored vector remains the target `pik`, while the true finite-population first-order inclusion probabilities are only approximately equal to that target. The five entry-point functions are organized by two dimensions, plus a dedicated dispatcher for balanced sampling: | | Without replacement | With replacement | |----------------|---------------------|------------------| | Equal prob. | `equal_prob_wor()` | `equal_prob_wr()` | | Unequal prob. | `unequal_prob_wor()`| `unequal_prob_wr()`| | Balanced | `balanced_wor()` | | ## Equal probability sampling ### Without replacement The simplest case: draw `n` units from a population of size `N` with equal probabilities. ```{r} library(sondage) s <- equal_prob_wor(N = 50, n = 5) s ``` The returned object stores the sampled unit indices and the inclusion probabilities: ```{r} s$sample s$pik[1:10] # all equal to n/N = 0.1 ``` Alternative methods are available. Systematic sampling provides implicit stratification based on the ordering of units: ```{r} s_sys <- equal_prob_wor(N = 50, n = 5, method = "systematic") s_sys ``` 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. ```{r} s_ber <- equal_prob_wor(N = 50, n = 5, method = "bernoulli") s_ber length(s_ber$sample) # may differ from 5 ``` ### With replacement ```{r} s_wr <- equal_prob_wr(N = 50, n = 5) s_wr ``` With-replacement designs track **hits** (how many times each unit was selected) instead of inclusion probabilities: ```{r} s_wr$hits # length-N vector of selection counts ``` ## Unequal probability sampling When units have different sizes (employees, revenue, area, ...), probability-proportional-to-size (PPS) designs give larger units higher selection probabilities, which reduces the variance of Horvitz-Thompson estimators [@horvitz1952]. ### Computing design quantities Start from a vector of positive size measures and compute inclusion probabilities (WOR) or expected hits (WR): ```{r} x <- c(10, 20, 5, 40, 25) n <- 3 # Without replacement: inclusion probabilities (capped at 1) pik <- inclusion_prob(x, n) pik sum(pik) # equals n # With replacement: expected hits (can exceed 1) eh <- expected_hits(x, n) eh sum(eh) # equals n ``` ### Without replacement Pass the inclusion probability vector to `unequal_prob_wor()`. The default method is Conditional Poisson Sampling [CPS, @chen1994], which has maximum entropy among fixed-size designs: ```{r} s <- unequal_prob_wor(pik, method = "cps") s s$sample inclusion_prob(s) # extract the stored design-defining target pik ``` Other fixed-size methods include `"brewer"` (draw-by-draw, order-invariant), `"sampford"` [@sampford1967], `"systematic"` (O(N), implicit stratification), `"sps"` [Sequential Poisson sampling, @ohlsson1998, supports PRN], and `"pareto"` [Pareto sampling, @rosen1997, 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()`. @chromy1979's method is the default: ```{r} s_wr <- unequal_prob_wr(eh, method = "chromy") s_wr s_wr$hits # realized selection counts ``` 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: ```{r} eh_big <- expected_hits(c(10, 20, 5, 40, 125), n) eh_big unequal_prob_wr(eh_big, method = "systematic_mr")$hits ``` ## Balanced sampling When auxiliary information is available at the design stage, balanced sampling ensures that Horvitz-Thompson estimates of auxiliary totals are (approximately) exact: $\sum_{k \in S} x_k / \pi_k \approx \sum_{k \in U} x_k$. This reduces variance for any study variable correlated with the balancing variables. The `balanced_wor()` function implements the cube method [@deville2004]: ```{r} pik <- inclusion_prob(c(10, 20, 5, 40, 25), n = 3) x <- matrix(c(10, 20, 5, 40, 25)) # auxiliary variable to balance on s_bal <- balanced_wor(pik, aux = x) s_bal ``` The cube method proceeds in two phases: a *flight phase* that moves probabilities toward 0 or 1 while maintaining all balancing constraints, and a *landing phase* that resolves the remaining undecided units by progressively relaxing constraints (starting from the last column of `aux`). The sample size constraint is always placed first and never relaxed. Users should order auxiliary variables by importance (most important first) in the `aux` matrix. For ill-conditioned or collinear auxiliary variables, `condition_aux = TRUE` applies weighted centering and scaling followed by rank pruning before running the cube algorithm. ### Stratified balanced sampling When strata are known, the stratified cube method [@chauvet2009] 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`): ```{r} 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 # Check per-stratum sample sizes table(strata[s_strat$sample]) ``` 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)` [@tripet2026]. 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: ```{r} 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]) ``` ### Spatially balanced sampling It's possible to select well-spread, spatially balanced samples using `lpm2` [local pivotal method 2, @grafstrom2012] and `scps` [spatially correlated Poisson sampling, @grafstrom2012scps]. ```{r} 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 s_scps$sample ``` 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 [@bondesson2008; @grafstrom2018]. Two surveys drawn with `prn = u` and `prn = 1 - u` are then negatively coordinated. @matei2023 spread on a measure of response burden instead of coordinates, so that highly burdened units rarely enter the same sample, and rarely enter both. ```{r} 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) ``` ## Design queries Once you have a design object, generics compute the quantities needed for variance estimation. ### Joint inclusion probabilities (WOR) ```{r} pik <- inclusion_prob(c(10, 20, 5, 40, 25), n = 3) s <- unequal_prob_wor(pik, method = "cps") pikl <- joint_inclusion_prob(s) pikl ``` 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) ```{r} eh <- expected_hits(c(10, 20, 5, 40, 25), n = 3) s_wr <- unequal_prob_wr(eh, method = "chromy") joint_expected_hits(s_wr) ``` ### 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: ```{r} 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) # Only the 10 x 10 submatrix for sampled units pikl_sub <- joint_inclusion_prob(s_large, sampled_only = TRUE) dim(pikl_sub) ``` 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`: ```{r} delta_sub <- sampling_cov(s_large, sampled_only = TRUE) dim(delta_sub) ``` ### 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): ```{r} sampling_cov(s) ``` With `weighted = TRUE`, it returns the Sen-Yates-Grundy check quantities $1 - \pi_i \pi_j / \pi_{ij}$ [@sen1953; @yatesandgrundy1953]. For well-behaved WOR designs, these should be non-positive off-diagonal: ```{r} sampling_cov(s, weighted = TRUE) ``` ## Batch sampling All five dispatchers accept an `nrep` parameter for drawing multiple independent samples in one call, useful for simulation studies. The result is always a design object (the same class as `nrep = 1`) but with `$sample` holding all replicates: ```{r} sim <- equal_prob_wor(N = 50, n = 5, nrep = 4) sim # design object dim(sim$sample) # 5 x 4 matrix: each column is one sample # Generics still work on batch objects inclusion_prob(sim) ``` For fixed-size designs `$sample` is a matrix (`n x nrep`). For random-size designs (Bernoulli, Poisson) it is a list of integer vectors, and realized sample sizes are available via `lengths(sim$sample)`. ## Design properties and statistical guarantees The table below summarizes what `sondage` computes for each method. Understanding these properties is essential for choosing the right variance estimator. ```{r, echo = FALSE} props <- data.frame( Method = c( "SRS WOR", "Systematic (EP)", "Bernoulli", "SRS WR", "CPS", "Brewer", "Sampford", "Systematic PPS", "Poisson", "SPS", "Pareto", "Chromy", "Multinomial", "Systematic PPS (MR)", "Cube", "LPM2", "SCPS" ), Type = c( rep("EP WOR", 3), "EP WR", rep("UP WOR", 7), rep("UP WR", 3), rep("Balanced", 3) ), `Fixed n` = c( "Yes", "Yes", "No", "Yes", "Yes", "Yes", "Yes", "Yes", "No", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes" ), `Exact design marginals?` = c( "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Target only*", "Target only*", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes" ), `Joint probs` = c( "Exact", "Exact (zeros)", "Exact", "Exact", "Exact", "HE approx", "Exact", "Exact (zeros)", "Exact", "HE approx", "HE approx", "Simulated", "Exact", "Exact (zeros)", "HE approx", "Not available", "Not available" ), check.names = FALSE ) knitr::kable(props, align = "llccl", caption = paste( "HE approx = high-entropy approximation (Brewer & Donadio, 2003).", "For WOR methods, design marginals are first-order inclusion probabilities pi_k,", "and for WR methods, they are expected hits E(N_k).", "Target only* = the stored design-defining pik is the target vector.", "For sps and pareto the true finite-population first-order probabilities", "are only approximately equal to that target, with the discrepancy", "vanishing asymptotically.", "Simulated = Monte Carlo via nsim parameter." ) ) ``` ### Notes on specific methods #### Equal-probability systematic sampling In `sondage`, `joint_inclusion_prob()` returns exact joint inclusion probabilities computed via C code, using the same algorithm as systematic PPS with uniform $\pi_k = n/N$. Some $\pi_{ij}$ are structurally zero (pairs of units that can never co-occur in the same systematic sample), so the SYG estimator is generally inapplicable. Consider replication or successive-differences variance estimators for systematic designs. #### SPS and Pareto (order sampling) Both are implemented via @rosen1997's order sampling framework. SPS uses ranking key $\xi_k = u_k / \pi_k$ [equivalent to @ohlsson1998's sequential threshold adjustment], and Pareto uses $\xi_k = [u_k/(1-u_k)] / [\pi_k/(1-\pi_k)]$. The $n$ units with the smallest $\xi_k$ are selected. Design objects store the supplied `pik` as the design-defining target vector, and `inclusion_prob(s)` returns that stored target. The true finite-population first-order inclusion probabilities are approximately (not exactly) equal to the target $\pi_k$ for finite populations, and the discrepancy vanishes as $N$ grows. Both support permanent random numbers (PRN) for sample coordination. #### Systematic PPS Joint inclusion probabilities are exact (computed via C code), but some $\pi_{ij}$ are structurally zero, for pairs of units that can never co-occur. `sampling_cov(s, weighted = TRUE)` returns `NA` for those entries and issues a warning, because the SYG estimator is undefined for such pairs. Consider successive-differences, Hartley-Rao, or replication methods instead. #### Chromy Pairwise expectations $E(n_i n_j)$ are estimated by Monte Carlo simulation in `joint_expected_hits()` (controlled by `nsim`, default 10 000). #### Systematic PPS with minimum replacement Each unit receives the whole part of its expected hits or one more, and the hits sum to $n$. Pairwise expectations from `joint_expected_hits()` are exact. A unit with at least one expected hit is selected in every sample, so its pairs with every unit of positive expected hits are positive. Two units below one expected hit can still be a structural zero, as in systematic PPS, and `sampling_cov(s, weighted = TRUE)` returns `NA` for those pairs. In a two-stage design, a primary unit hit twice needs two independent second-stage samples, and the with-replacement variance approximation then treats each hit as a separate primary unit. #### Cube (balanced sampling) For the cube method [@deville2004], `joint_inclusion_prob()` uses the high-entropy approximation. It depends only on `pik`, so it does not reflect strata or balancing constraints. A stratified design can have pairs that never occur together, and the approximation reports a positive probability for them. For a stratified or strongly constrained design, prefer Monte Carlo estimates of the joint probabilities or a variance method suited to balanced sampling. The stratified variant [@chauvet2009] preserves exact within-stratum sample sizes. Order auxiliary variables by importance (most important first) since the landing phase relaxes constraints from last to first. Inequality bounds provide controlled selection, but tight bounds can make the high-entropy approximation less accurate. #### LPM2 and SCPS (spatially balanced sampling) Both methods preserve exact first-order inclusion probabilities and spread the sample using the variables supplied in `spread`. Their low-entropy spatial dependence is deliberate, so no joint inclusion probability or sampling covariance approximation is exposed. Use a local-neighborhood variance estimator instead. ### 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 [@sen1953; @yatesandgrundy1953] 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 [@brewer2003, 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 @hansen1943 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 @chromy2009 which unifies WOR and WR variance estimation via `sampling_cov()`. ## Extending sondage `register_method()` lets you plug custom unequal-probability sampling algorithms into the existing dispatchers and generics. Registered methods work with `unequal_prob_wor()`, `joint_inclusion_prob()`, `sampling_cov()`, batch mode, and all downstream tooling. See `vignette("custom-methods")` for worked examples. ## References