A Tidy Grammar for Survey Sampling • samplyr

A tidy grammar for survey sampling in R. samplyr provides a minimal set of composable verbs for stratified, clustered, multi-stage, and multi-phase sampling designs with PPS methods, sample coordination, and panel rotation.

Ecosystem Positioning

samplyr sits in the middle of a three-layer workflow:

Use svyplan for planning (n, precision, power, allocation, budget).
Use samplyr to specify designs, draw samples, and carry design metadata through fieldwork.
Use survey/srvyr after data collection for estimation and inference.

These packages are complementary, not competing. Handoffs are explicit (for example, svyplan outputs can feed draw(n = ...), and samplyr outputs convert via as_svydesign(), as_svrepdesign(), and as_survey_design()).

Learn samplyr

Follow two articles: Get started for the grammar and a complete three-stage fieldwork example, then Analyze the sample for the handoff to survey and srvyr.

The book Survey Sampling Design with R covers planning with svyplan, coordination, rotating panels and validation. The selection_methods reference compares the methods, and the variance_estimation reference gives how each one’s variance is estimated. Serialization and design semantics are references to consult as needed.

Why samplyr?

samplyr is built around a simple idea: sampling code should read like its English description.

library(samplyr)
data(bfa_eas)

# "Stratify by region, proportionally allocate 500 samples, execute"
overview_sample <- sampling_design() |>
  stratify_by(region, alloc = "proportional") |>
  draw(n = 500) |>
  execute(bfa_eas, seed = 1)

overview_sample |>
  as.data.frame() |>
  dplyr::select(ea_id, region, households, .weight) |>
  head(4)
#>   ea_id            region households  .weight
#> 1  9648 Boucle du Mouhoun         35 89.44643
#> 2 11547 Boucle du Mouhoun          7 89.44643
#> 3 41824 Boucle du Mouhoun         10 89.44643
#> 4 11012 Boucle du Mouhoun         71 89.44643

Consider a real survey design from Lohr (2022, Example 7.1), based on a 1991 study of bed net use in rural Gambia (D’Alessandro et al., 1994):

Malaria morbidity can be reduced by using bed nets impregnated with insecticide, but this is only effective if the bed nets are in widespread use. In 1991, a nationwide survey was designed to estimate the prevalence of bed net use in rural areas of the Gambia (D’Alessandro et al., 1994).

The sampling frame consisted of all rural villages of fewer than 3,000 people. The villages were stratified by three geographic regions (eastern, central, and western) and by whether the village had a public health clinic (PHC) or not. In each region five districts were chosen with probability proportional to the district population. In each district four villages were chosen, again with probability proportional to census population: two PHC villages and two non-PHC villages. Finally, six compounds were chosen more or less randomly from each village.

In samplyr, this three-stage stratified cluster design translates directly into code:

design <- sampling_design(title = "Gambia bed nets") |>
  add_stage() |>
    stratify_by(region) |>
    cluster_by(district) |>
    draw(n = 5, method = "pps_brewer", mos = district_population) |>
  add_stage() |>
    stratify_by(phc) |>
    cluster_by(village) |>
    draw(n = 2, method = "pps_brewer", mos = village_population) |>
  add_stage() |>
    draw(n = 6)
design
#> ── Sampling Design: Gambia bed nets ─────────────────────────────
#> 
#> ℹ 3 stages
#> 
#> ── Stage 1 ──────────────────────────────────────────────────────
#> • Strata: region
#> • Cluster: district
#> • Draw: n = 5 (per stratum), method = pps_brewer, mos = district_population
#> 
#> ── Stage 2 ──────────────────────────────────────────────────────
#> • Strata: phc
#> • Cluster: village
#> • Draw: n = 2 (per stratum, per district), method = pps_brewer, mos = village_population
#> 
#> ── Stage 3 ──────────────────────────────────────────────────────
#> • Draw: n = 6 (per village), method = srswor

The two MOS names distinguish district totals from village totals. This is a design declaration. Execution also needs compatible registers and compound listings. For a complete runnable example using bundled data, see A three-stage household sample in vignette("introduction"): communes, EAs and synthetic household listings, followed by observations and survey export.

Lohr, S. L. (2022). Sampling: Design and Analysis (3rd ed.). CRC Press.

Installation

# sondage, the sampling algorithms backend, is on CRAN
install.packages("sondage")

# Install svyplan (sample size, precision, power, and stratification)
pak::pkg_install("gitlab::dickoa/svyplan")

# Install samplyr
pak::pkg_install("gitlab::dickoa/samplyr")

The Grammar

samplyr uses 5 verbs and 1 modifier:

Function	Purpose
`sampling_design()`	Create a new sampling design
`stratify_by()`	Define stratification and allocation
`cluster_by()`	Define cluster/PSU variable
`draw()`	Specify sample size and method
`execute()`	Run the design on a frame
`add_stage()`	Delimit stages in multi-stage designs

Frame-Independent Design

stratify_by() and cluster_by() take bare column names. The design is stored as a specification and resolved only when a frame is available (validate_frame(), execute(), as_svydesign()), so design specification stays separate from execution.

design <- sampling_design() |>
  stratify_by(region, alloc = "proportional") |>
  cluster_by(ea_id) |>
  draw(n = 300)

sample <- execute(design, bfa_eas, seed = 2)
sample
#> # A tbl_sample: 300 × 18
#> # Sampling:     1 stage | 300 units
#> # Weights:      148.57 [146.5, 151.22]
#>    ea_id region           province commune urban_rural population
#>  * <int> <fct>            <fct>    <fct>   <fct>            <int>
#>  1 29513 Boucle du Mouho… Bale     Bana    Rural              181
#>  2 29527 Boucle du Mouho… Bale     Bana    Rural              193
#>  3 36703 Boucle du Mouho… Bale     Fara    Rural              975
#>  4  8455 Boucle du Mouho… Bale     Ouri    Rural               57
#>  5  8580 Boucle du Mouho… Bale     Pâ      Rural              171
#>  6 11739 Boucle du Mouho… Bale     Yaho    Rural              509
#>  7 11746 Boucle du Mouho… Bale     Yaho    Rural               98
#>  8  6291 Boucle du Mouho… Banwa    Kouka   Rural              563
#>  9 34031 Boucle du Mouho… Banwa    Sami    Rural               41
#> 10 34058 Boucle du Mouho… Banwa    Sami    Rural              516
#> # ℹ 290 more rows
#> # ℹ 12 more variables: households <int>, area_km2 <dbl>,
#> #   pop_density <dbl>, longitude <dbl>, latitude <dbl>,
#> #   remoteness <fct>, fieldwork_cost <int>, .weight <dbl>,
#> #   .sample_id <int>, .stage <int>, .weight_1 <dbl>,
#> #   .fpc_1 <dbl>

Quick Start

library(samplyr)
data(bfa_eas)

# Simple random sample
srs_smpl <- sampling_design() |>
  draw(n = 100) |>
  execute(bfa_eas, seed = 321)

srs_smpl
#> # A tbl_sample: 100 × 18
#> # Sampling:     1 stage | 100/44,570 units
#> # Weights:      445.7 [445.7, 445.7]
#>    ea_id region           province commune urban_rural population
#>  * <int> <fct>            <fct>    <fct>   <fct>            <int>
#>  1  5168 Centre-Nord      Sanmate… Kaya    Urban              794
#>  2 11281 Centre-Sud       Nahouri  Tiébélé Rural              123
#>  3 30149 Centre-Ouest     Sanguie  Dassa   Rural              813
#>  4 13590 Est              Gnagna   Koala   Rural              353
#>  5 26045 Boucle du Mouho… Mouhoun  Dédoug… Rural              378
#>  6 44140 Centre-Est       Koulpel… Soudou… Rural              258
#>  7 30931 Est              Komandj… Gayéri  Rural              150
#>  8 10745 Sahel            Yagha    Tankou… Rural              312
#>  9 13990 Hauts-Bassins    Kenedou… Kourou… Rural              406
#> 10 27757 Est              Tapoa    Partia… Rural              571
#> # ℹ 90 more rows
#> # ℹ 12 more variables: households <int>, area_km2 <dbl>,
#> #   pop_density <dbl>, longitude <dbl>, latitude <dbl>,
#> #   remoteness <fct>, fieldwork_cost <int>, .weight <dbl>,
#> #   .sample_id <int>, .stage <int>, .weight_1 <dbl>,
#> #   .fpc_1 <int>

# Stratified proportional allocation
strata_smpl <- sampling_design() |>
  stratify_by(region, alloc = "proportional") |>
  draw(n = 300) |>
  execute(bfa_eas, seed = 12)

strata_smpl
#> # A tbl_sample: 300 × 18
#> # Sampling:     1 stage | 300/44,570 units
#> # Weights:      148.57 [146.5, 151.22]
#>    ea_id region           province commune urban_rural population
#>  * <int> <fct>            <fct>    <fct>   <fct>            <int>
#>  1 31987 Boucle du Mouho… Sourou   Kiemba… Rural               24
#>  2  4968 Boucle du Mouho… Sourou   Kassoum Rural              303
#>  3  4958 Boucle du Mouho… Sourou   Kassoum Rural              142
#>  4 23903 Boucle du Mouho… Banwa    Solenzo Rural               94
#>  5 25996 Boucle du Mouho… Mouhoun  Dédoug… Rural              538
#>  6 21501 Boucle du Mouho… Kossi    Doumba… Rural              247
#>  7 43794 Boucle du Mouho… Mouhoun  Safané  Rural               87
#>  8  9724 Boucle du Mouho… Banwa    Sanaba  Rural               51
#>  9 10979 Boucle du Mouho… Banwa    Tansila Rural               71
#> 10  8893 Boucle du Mouho… Bale     Pompoï  Rural              522
#> # ℹ 290 more rows
#> # ℹ 12 more variables: households <int>, area_km2 <dbl>,
#> #   pop_density <dbl>, longitude <dbl>, latitude <dbl>,
#> #   remoteness <fct>, fieldwork_cost <int>, .weight <dbl>,
#> #   .sample_id <int>, .stage <int>, .weight_1 <dbl>,
#> #   .fpc_1 <dbl>

# PPS cluster sampling
cluster_smpl <- sampling_design() |>
  cluster_by(ea_id) |>
  draw(n = 50, method = "pps_brewer", mos = households) |>
  execute(bfa_eas, seed = 123)

cluster_smpl
#> # A tbl_sample: 50 × 19
#> # Sampling:     1 stage | 50 units
#> # Weights:      779.3 [188.64, 2999.41]
#>    ea_id region           province commune urban_rural population
#>  * <int> <fct>            <fct>    <fct>   <fct>            <int>
#>  1 23944 Boucle du Mouho… Banwa    Solenzo Rural              758
#>  2 33134 Boucle du Mouho… Kossi    Nouna   Rural             1116
#>  3 10578 Boucle du Mouho… Kossi    Sônô    Rural              326
#>  4 23289 Cascades         Comoe    Ouô     Rural              534
#>  5 15240 Centre           Kadiogo  Ouagad… Urban              614
#>  6 15327 Centre           Kadiogo  Ouagad… Urban             1061
#>  7 15429 Centre           Kadiogo  Ouagad… Urban              667
#>  8 15550 Centre           Kadiogo  Ouagad… Urban              967
#>  9 17098 Centre           Kadiogo  Ouagad… Urban              844
#> 10 17487 Centre           Kadiogo  Ouagad… Urban              938
#> # ℹ 40 more rows
#> # ℹ 13 more variables: households <int>, area_km2 <dbl>,
#> #   pop_density <dbl>, longitude <dbl>, latitude <dbl>,
#> #   remoteness <fct>, fieldwork_cost <int>, .weight <dbl>,
#> #   .sample_id <int>, .stage <int>, .weight_1 <dbl>,
#> #   .fpc_1 <int>, .certainty_1 <lgl>

Multi-Stage Sampling

Use add_stage() to define multi-stage designs. This example selects districts with PPS, then samples EAs within each:

library(dplyr, warn.conflicts = FALSE)
data(zwe_eas)

# Add district-level measure of size
zwe_frame <- zwe_eas |>
  mutate(district_hh = sum(households), .by = district)

# Two-stage design: 10 districts, 5 EAs per district
sample <- sampling_design() |>
  add_stage(label = "Districts") |>
    cluster_by(district) |>
    draw(n = 10, method = "pps_brewer", mos = district_hh) |>
  add_stage(label = "EAs") |>
    draw(n = 5) |>
  execute(zwe_frame, seed = 12345)
#> Stage 1: 1 unit is selected with certainty.
#> ℹ A unit whose share of its pool's size is above 1 / n has
#>   chance 1, and the rest of the take is spread over the other
#>   units.
#> ℹ Set `certainty_size` or `certainty_prop` in `draw()` to state
#>   the rule in the design, or check before drawing with
#>   `frame_summary()`.

sample
#> # A tbl_sample: 50 × 21
#> # Sampling:     2 stages | 50/107,250 units
#> # Weights:      2514.13 [784.4, 4121.21]
#>    ea_id province   district    ward_pcode urban_rural population
#>  * <int> <fct>      <fct>       <chr>      <fct>            <int>
#>  1 88821 Harare     Harare      ZW192102   Urban              197
#>  2 87420 Harare     Harare      ZW192116   Urban             1378
#>  3 86941 Harare     Harare      ZW192108   Urban              103
#>  4 88444 Harare     Harare      ZW192129   Urban              538
#>  5 34696 Harare     Harare      ZW192120   Urban              452
#>  6 14864 Manicaland Chimanimani ZW110201   Rural              111
#>  7 59668 Manicaland Chimanimani ZW110204   Rural              142
#>  8 60166 Manicaland Chimanimani ZW110203   Rural              114
#>  9 60082 Manicaland Chimanimani ZW110216   Rural               73
#> 10 58991 Manicaland Chimanimani ZW110219   Rural               73
#> # ℹ 40 more rows
#> # ℹ 15 more variables: households <int>, buildings <int>,
#> #   women_15_49 <int>, men_15_49 <int>, children_under5 <int>,
#> #   area_km2 <dbl>, district_hh <int>, .weight <dbl>,
#> #   .sample_id <int>, .stage <int>, .weight_2 <dbl>,
#> #   .fpc_2 <dbl>, .weight_1 <dbl>, .fpc_1 <int>,
#> #   .certainty_1 <lgl>

Operational Sampling

Execute stages separately when fieldwork happens between stages. The example follows a common household survey protocol, with EAs selected by systematic PPS on a frame sorted by province and district, and a fixed take of households selected systematically from the new listing.

design <- sampling_design() |>
  add_stage(label = "EA") |>
    stratify_by(urban_rural) |>
    cluster_by(ea_id) |>
    draw(n = 10, method = "pps_systematic", mos = households,
         control = c(province, district)) |>
  add_stage(label = "HH") |>
    draw(n = 5, method = "systematic")

# Execute stage 1 only
selected_eas <- execute(design, zwe_eas, stages = 1, seed = 1)
selected_eas
#> # A tbl_sample: 20 × 18
#> # Stages:       1/2
#> # Sampling:     1 stage | 20 units
#> # Weights:      4121.95 [168.15, 13526.55]
#>     ea_id province     district ward_pcode urban_rural population
#>  *  <int> <fct>        <fct>    <chr>      <fct>            <int>
#>  1  20129 Bulawayo     Bulawayo ZW102114   Urban              449
#>  2  85826 Harare       Chitung… ZW192204   Urban             4120
#>  3  85720 Harare       Harare   ZW192144   Urban             1012
#>  4  87853 Harare       Harare   ZW192145   Urban              656
#>  5  33471 Harare       Harare … ZW192401   Urban              140
#>  6  57994 Manicaland   Chimani… ZW110220   Rural              173
#>  7  71844 Manicaland   Mutare … ZW112109   Urban              375
#>  8  31050 Manicaland   Mutasa   ZW110619   Rural               95
#>  9  95223 Mashonaland… Mazowe   ZW120412   Rural              191
#> 10  83444 Mashonaland… Goromon… ZW130216   Urban              321
#> 11  84047 Mashonaland… Maronde… ZW130403   Rural              138
#> 12  81288 Mashonaland… Chegutu  ZW140125   Rural              111
#> 13  35868 Mashonaland… Chinhoyi ZW142113   Urban             1785
#> 14  33035 Mashonaland… Zvimba   ZW140635   Rural              286
#> 15  27459 Masvingo     Chiredz… ZW182208   Urban              358
#> 16  56591 Masvingo     Gutu     ZW180408   Rural              106
#> 17  31399 Matabelelan… Hwange   ZW150309   Rural               67
#> 18 102583 Matabelelan… Insiza   ZW160513   Rural              103
#> 19  69632 Midlands     Gokwe S… ZW172903   Urban              717
#> 20  10452 Midlands     Kwekwe   ZW170531   Rural              120
#> # ℹ 12 more variables: households <int>, buildings <int>,
#> #   women_15_49 <int>, men_15_49 <int>, children_under5 <int>,
#> #   area_km2 <dbl>, .weight <dbl>, .sample_id <int>,
#> #   .stage <int>, .weight_1 <dbl>, .fpc_1 <int>,
#> #   .certainty_1 <lgl>

# ... fieldwork ...
# list all households in selected EAs
selected_eas_list <- selected_eas |>
  slice(rep(seq_len(n()), households)) |>
  mutate(hh_id = row_number())

# Execute stage 2
final_sample <- selected_eas |>
  execute(selected_eas_list, seed = 2)
final_sample
#> # A tbl_sample: 100 × 21
#> # Sampling:     2 stages | 100 units
#> # Weights:      38147.25 [35714.84, 40579.66]
#>    ea_id province district    ward_pcode urban_rural population
#>  * <int> <fct>    <fct>       <chr>      <fct>            <int>
#>  1 20129 Bulawayo Bulawayo    ZW102114   Urban              449
#>  2 20129 Bulawayo Bulawayo    ZW102114   Urban              449
#>  3 20129 Bulawayo Bulawayo    ZW102114   Urban              449
#>  4 20129 Bulawayo Bulawayo    ZW102114   Urban              449
#>  5 20129 Bulawayo Bulawayo    ZW102114   Urban              449
#>  6 85826 Harare   Chitungwiza ZW192204   Urban             4120
#>  7 85826 Harare   Chitungwiza ZW192204   Urban             4120
#>  8 85826 Harare   Chitungwiza ZW192204   Urban             4120
#>  9 85826 Harare   Chitungwiza ZW192204   Urban             4120
#> 10 85826 Harare   Chitungwiza ZW192204   Urban             4120
#> # ℹ 90 more rows
#> # ℹ 15 more variables: households <int>, buildings <int>,
#> #   women_15_49 <int>, men_15_49 <int>, children_under5 <int>,
#> #   area_km2 <dbl>, hh_id <int>, .weight <dbl>,
#> #   .sample_id <int>, .stage <int>, .weight_2 <dbl>,
#> #   .fpc_2 <int>, .weight_1 <dbl>, .fpc_1 <int>,
#> #   .certainty_1 <lgl>

selected_eas still contains design plus the realized stage-1 selection. Using it as the first argument continues that same multi-stage design. The expanded listing is only the candidate frame for stage 2. Starting from a new design with a prior tbl_sample as its frame instead denotes a new sampling phase.

If tidyr::uncount() or another operation drops the listing’s tbl_sample class, keep using that plain object as the second argument above. execute() recognizes its inherited sampling attributes and generated columns and refuses to run it as a fresh frame from design, which would silently rerun stage 1. Passing the intact selected_eas back as a frame of design also warns: that call is supported as a new sampling phase, but it is not stage continuation.

Selection Methods

Sixteen methods ship in three families:

Equal probability: srswor (the default), srswr, systematic, bernoulli
PPS, all requiring a measure of size: pps_brewer, pps_systematic, pps_cps, pps_sampford, pps_poisson, pps_sps, pps_pareto, pps_multinomial, pps_chromy
Balanced: cube, plus the spatial lpm2 and scps

?selection_methods is the reference: which take n or frac, which need mos, aux or spread, which have a random rather than fixed sample size, which accept permanent random numbers, and the paper behind each. Custom methods registered with sondage::register_method() are used the same way.

Allocation Methods

When stratifying, n in draw() has one of two meanings. Without an allocation method it is the size taken in each stratum. With one, it is the total for the stage, shared across the strata by the method.

Method	Description
(none)	`n` is taken in each stratum
`equal`	Splits the total `n` as evenly as the strata allow
`proportional`	Proportional to stratum size
`neyman`	Minimize variance (requires `variance`)
`optimal`	Minimize cost-variance (requires `variance` and `cost`)
`power`	Compromise allocation (requires `cv` and `importance`)

The same n = 8 gives 24 units without a method and 8 with one.

frame <- data.frame(region = rep(c("A", "B", "C"), each = 10), id = 1:30)

per_stratum <- sampling_design() |>
  stratify_by(region) |>
  draw(n = 8) |>
  execute(frame, seed = 1)
table(per_stratum$region)
#> 
#> A B C 
#> 8 8 8

total <- sampling_design() |>
  stratify_by(region, alloc = "equal") |>
  draw(n = 8) |>
  execute(frame, seed = 1)
table(total$region)
#> 
#> A B C 
#> 3 3 2

Shares are rounded by largest remainder, which keeps the total. Equal remainders go to the larger target first, then in label order, so the three equal shares of 8 become 3, 3 and 2. validate_frame() states the split before anything is drawn, and the tied column of frame_summary(detail = "pool") marks the strata a tie decided.

At a later stage, both readings apply within each selected parent. A total at stage 2 is a total per selected district, not one for the whole survey. The printed design states the scope of every n.

sampling_design() |>
  add_stage(label = "District") |>
    stratify_by(region) |>
    cluster_by(district) |>
    draw(n = 5) |>
  add_stage(label = "Village") |>
    stratify_by(urban, alloc = "equal") |>
    cluster_by(village) |>
    draw(n = 4)
#> ── Sampling Design ──────────────────────────────────────────────
#> 
#> ℹ 2 stages
#> 
#> ── Stage 1: District ────────────────────────────────────────────
#> • Strata: region
#> • Cluster: district
#> • Draw: n = 5 (per stratum), method = srswor
#> 
#> ── Stage 2: Village ─────────────────────────────────────────────
#> • Strata: urban (equal)
#> • Cluster: village
#> • Draw: n = 4 (total, per district), method = srswor

Sample Size Bounds

Use min_n and max_n in draw() to constrain stratum sample sizes when using allocation methods:

data(bfa_eas_variance)

# Ensure at least 2 per stratum (minimum for variance estimation)
sampling_design() |>
  stratify_by(region, alloc = "neyman", variance = bfa_eas_variance) |>
  draw(n = 300, min_n = 2) |>
  execute(bfa_eas, seed = 321)
#> # A tbl_sample: 300 × 18
#> # Sampling:     1 stage | 300/44,570 units
#> # Weights:      148.57 [142.32, 162]
#>    ea_id region           province commune urban_rural population
#>  * <int> <fct>            <fct>    <fct>   <fct>            <int>
#>  1 11643 Boucle du Mouho… Nayala   Yaba    Rural              396
#>  2 36763 Boucle du Mouho… Bale     Fara    Rural               79
#>  3  9648 Boucle du Mouho… Banwa    Sanaba  Rural              279
#>  4 10555 Boucle du Mouho… Kossi    Sônô    Rural              164
#>  5 34926 Boucle du Mouho… Sourou   Tougan  Rural             1239
#>  6 44501 Boucle du Mouho… Mouhoun  Tchéri… Rural              221
#>  7  9702 Boucle du Mouho… Banwa    Sanaba  Rural              260
#>  8 26045 Boucle du Mouho… Mouhoun  Dédoug… Rural              378
#>  9  8605 Boucle du Mouho… Bale     Pâ      Rural              283
#> 10 21078 Boucle du Mouho… Kossi    Bouras… Rural               93
#> # ℹ 290 more rows
#> # ℹ 12 more variables: households <int>, area_km2 <dbl>,
#> #   pop_density <dbl>, longitude <dbl>, latitude <dbl>,
#> #   remoteness <fct>, fieldwork_cost <int>, .weight <dbl>,
#> #   .sample_id <int>, .stage <int>, .weight_1 <dbl>,
#> #   .fpc_1 <dbl>

Custom Allocation

For custom stratum-specific sizes or rates, pass a data frame to n or frac in draw():

sizes_df <- data.frame(
  region = c("A", "B", "C"),
  n = c(2, 3, 4)
)
## sizes <- c(A = 2, B = 3, C = 4) can also be used with one stratum variable

sample <- sampling_design() |>
  stratify_by(region) |>
  draw(n = sizes_df) |>
  execute(frame, seed = 42)
table(sample$region)
#> 
#> A B C 
#> 2 3 4

# Neyman allocation
data(bfa_eas_variance)

sample <- sampling_design() |>
  stratify_by(region, alloc = "neyman", variance = bfa_eas_variance) |>
  draw(n = 300) |>
  execute(bfa_eas, seed = 4321)

sample
#> # A tbl_sample: 300 × 18
#> # Sampling:     1 stage | 300/44,570 units
#> # Weights:      148.57 [142.32, 162]
#>    ea_id region           province commune urban_rural population
#>  * <int> <fct>            <fct>    <fct>   <fct>            <int>
#>  1 21517 Boucle du Mouho… Kossi    Doumba… Rural              354
#>  2 33112 Boucle du Mouho… Kossi    Nouna   Rural               28
#>  3 25963 Boucle du Mouho… Mouhoun  Dédoug… Rural              416
#>  4  6354 Boucle du Mouho… Banwa    Kouka   Rural              296
#>  5 11787 Boucle du Mouho… Bale     Bagassi Rural              874
#>  6  8233 Boucle du Mouho… Mouhoun  Ouarko… Rural              183
#>  7  9925 Boucle du Mouho… Bale     Siby    Rural              204
#>  8 34812 Boucle du Mouho… Sourou   Tougan  Rural             1072
#>  9 11441 Boucle du Mouho… Sourou   Toéni   Rural               64
#> 10  8883 Boucle du Mouho… Bale     Pompoï  Rural              192
#> # ℹ 290 more rows
#> # ℹ 12 more variables: households <int>, area_km2 <dbl>,
#> #   pop_density <dbl>, longitude <dbl>, latitude <dbl>,
#> #   remoteness <fct>, fieldwork_cost <int>, .weight <dbl>,
#> #   .sample_id <int>, .stage <int>, .weight_1 <dbl>,
#> #   .fpc_1 <dbl>

Surveys that publish an estimate for every region often take a compromise between proportional and equal allocation (Bankier 1988). Power allocation with a common cv and the stratum size as importance gives each region a take proportional to the square root of its size.

ea_counts <- count(bfa_eas, region, name = "importance")

compromise <- sampling_design() |>
  stratify_by(region, alloc = "power",
              cv = data.frame(region = ea_counts$region, cv = 1),
              importance = ea_counts, power = 0.5) |>
  draw(n = 300) |>
  execute(bfa_eas, seed = 5)
table(compromise$region)
#> 
#> Boucle du Mouhoun          Cascades            Centre 
#>                28                20                25 
#>        Centre-Est       Centre-Nord      Centre-Ouest 
#>                22                23                24 
#>        Centre-Sud               Est     Hauts-Bassins 
#>                16                30                28 
#>              Nord   Plateau-Central             Sahel 
#>                22                16                26 
#>         Sud-Ouest 
#>                20

Beyond the Basics

The verbs above compose across many common design families. These capabilities extend the same grammar. The linked references describe their inputs and limits:

Capability	How	Where
Balanced sampling	`draw(method = "cube", aux = ...)`, with `bound()` for hard count constraints	`?draw`, `?bound`
Spatially balanced	`draw(method = "lpm2" or "scps", spread = c(lon, lat))`	`?selection_methods`
Sample coordination	`draw(prn = ...)` with permanent random numbers, for overlap across waves	`?draw`, Coordinate samples with permanent random numbers
Custom methods	`sondage::register_method()`, then `pps_<name>` or `balanced_<name>`	`?draw`, `sondage::register_method()`
Panel rotation	`execute(panels = 4)`, or `panel_stage =` to rotate units inside retained parents	`?panel_assignment`, `?rotation_program`, Plan and run a rotating panel
Replicated draws	`execute(reps = 5)`	`?execute`
Two-phase	pipe a `tbl_sample` into a new design’s `execute()`	`vignette("survey-analysis")`
Indirect sampling	`share_weights()`, to estimate for a population linked to the one that was sampled	`vignette("design-semantics")`
Overlapping frames	`stack_frames()`, for two registers of one population	`vignette("survey-analysis")`

The variance_estimation reference (?variance_estimation) compares first-order and joint probability quality and variance estimators, with the stage and phase restrictions.

# Balanced on auxiliary totals, PPS on size
sampling_design() |>
  stratify_by(region, alloc = "proportional") |>
  draw(n = 300, method = "cube", mos = households,
       aux = c(population, area_km2)) |>
  execute(bfa_eas, seed = 24)

# Four rotation groups, five independent replicates
execute(design, bfa_eas, seed = 1, panels = 4)
execute(design, bfa_eas, seed = 42, reps = 5)

# Areas stay in the survey while the households inside them rotate
execute(two_stage_design, frame, seed = 1, panels = 4, panel_stage = 2)

Survey Export

Convert to survey or srvyr for estimation. Here the Neyman sample above estimates the mean number of households per EA:

svy <- as_svydesign(sample)
survey::svymean(~households, svy)
#>              mean     SE
#> households 64.165 2.8661

srvyr::as_survey_design(sample) |>
  dplyr::summarise(mean_households = srvyr::survey_mean(households))
#> # A tibble: 1 × 2
#>   mean_households mean_households_se
#>             <dbl>              <dbl>
#> 1            64.2               2.87

For fixed-size PPS WOR stages, as_svydesign() uses the Brewer variance approximation by default. Other selection families receive their own treatment. joint_expectation() supplies exact or approximate pairwise quantities where supported. An exact matrix does not guarantee a stable variance estimate or accurate interval coverage. See vignette("survey-analysis") for domain estimation, replicate weights, and the method-by-method breakdown.

Both export verbs also take a stack_frames() collection, compositing overlapping frames through survey::multiframe() or through a combined replicate design that varies one frame at a time.

Diagnostics

summary(strata_smpl)
#> ── Sample Summary ───────────────────────────────────────────────
#> 
#> ℹ n = 300 of 44,570 | stages = 1/1 | seed = 12
#> 
#> ── Stage 1 ──────────────────────────────────────────────────────
#> • srswor, by region (proportional)
#> • 13 strata: N_h 1,612-5,505, n_h 11-37, f_h 0.0066-0.0068
#> 
#> ── Weights ──────────────────────────────────────────────────────
#> • Mean 148.57 [146.5, 151.22] | CV 0.01 | Kish DEFF 1 | n_eff 300

By default, each execution records a summary frame digest: a compact record of the selection pools and chances the design resolved. frame_summary() returns it as tibbles without needing the frame. Use frame_digest = "none" for the lowest execution overhead, or frame_digest = "full" when downstream work needs exact per-unit chances for unequal-probability element stages. The companion package samplens, in development, draws the digest as a visual sampling card. See vignette("introduction") for examples.

execute() also reports when the frame could not deliver what the design asked for. A pool holding fewer units than the stage requested is selected whole, which makes the design non-self-weighting, so this is a warning rather than a silent adjustment. The capped column of frame_summary(sample, detail = "pool") names the affected pools. ?execution_conditions lists every such condition and the payload it carries.

Checking a Design Before Drawing

A design can be checked against its frame before any unit is drawn, and none of the checks below uses random numbers. A measure of size belongs to the sampling unit, so it may repeat on every row of a cluster but must be the same on each. At the final stage the frame rows are the units and need no cluster_by().

villages <- data.frame(
  village = rep(c("v1", "v2", "v3", "v4"), c(6, 8, 5, 9)),
  household = 1:28
)
villages$size <- ave(villages$household, villages$village, FUN = length)

plan <- sampling_design() |>
  add_stage(label = "Village") |>
    cluster_by(village) |>
    draw(n = 2, method = "pps_brewer", mos = size) |>
  add_stage(label = "Household") |>
    draw(n = 6)

validate_frame(plan, villages)
#> Stage 2: the take exceeds the pool population in 1 of 4
#> candidate pools.
#> ℹ Short pools: "v3".
#> ℹ `execute()` takes every unit of a short pool and warns.
#>   Inspect with `frame_summary(design, frame, detail = "pool")`.
frame_summary(plan, villages)[c("stage", "n_pools", "n_capped", "n_target")]
#> # A tibble: 2 × 4
#>   stage n_pools n_capped n_target
#>   <int>   <int>    <int>    <dbl>
#> 1     1       1        0        2
#> 2     2       4        1       12

The pool_key column of frame_summary(plan, villages, detail = "pool") names the short pools as execute() will. By default execute() takes every unit of a short pool and warns. With on_capped = "error" in the stage’s draw(), it refuses before anything is drawn.

strict <- sampling_design() |>
  add_stage(label = "Village") |>
    cluster_by(village) |>
    draw(n = 2, method = "pps_brewer", mos = size) |>
  add_stage(label = "Household") |>
    draw(n = 6, on_capped = "error")

execute(strict, villages, seed = 1)
#> Error in `execute()`:
#> ! Stage 2: the take exceeds the pool population in 1 of 4
#>   candidate pools.
#> ℹ Short pools: "v3".
#> ℹ The stage sets `on_capped = "error"`, so nothing was drawn.
#> ℹ Lower the take, merge the short pools in the frame, or set
#>   `on_capped = "warn"` to take every unit of a short pool.

A requirement of the protocol, such as two villages and five households in each, can be written as a design with the default method and compared with the design to field. compare_designs() matches pools by the frame rows they hold and compares each pool’s take. The selection methods and measures of size are listed side by side but do not enter the verdict. With level = "probabilities" it compares the inclusion probability of every unit instead.

requirement <- sampling_design() |>
  add_stage() |>
    cluster_by(village) |>
    draw(n = 2) |>
  add_stage() |>
    draw(n = 5)

compare_designs(requirement, plan, villages, key = household)
#> ── Design comparison ────────────────────────────────────────────
#> Allocations differ on this frame.
#> 
#> ✔ Stage 1 units: agree. 4 units.
#> ✔ Stage 1 pools: agree. 1 pool.
#> ✔ Stage 1 targets: agree. Every pool takes the same. x takes 2 units in all, y 2.
#> ✔ Stage 2 units: agree. 28 units.
#> ✔ Stage 2 pools: agree. 4 pools.
#> ✖ Stage 2 targets: differ. 4 of 4 pools differ ("v1", "v2", "v3" and 1 more).
#> 
#> Methods and measures of size, outside the verdict:
#> ℹ Stage 1 method: srswor in x, pps_brewer in y.
#> ℹ Stage 1 measure of size: none in x, size in y.
#> ✔ Stage 2 method: srswor in both.

The plan takes six households in each village where the requirement asks for five, and the comparison names the villages.

Validation

For statistical validation on synthetic populations with known truths, see Validation on synthetic populations in the book. It combines deterministic invariants (weights, certainty, stage compounding) with Monte Carlo checks of bias, standard-error calibration, and coverage.

Included Datasets

Derived and synthetic sampling frames for learning and testing:

Dataset	Description	Rows
`bfa_eas`	Household budget and living-standards EA frame (Burkina Faso)	44,570
`zwe_eas`	Demographic, health, and child-indicator EA frame (Zimbabwe)	107,250
`ken_enterprises`	Establishment survey frame (Kenya)	17,004

Plus auxiliary data: bfa_eas_variance, bfa_eas_cost

Comparison with SAS and SPSS

Sampford’s method is the PPS method without replacement that the three tools share, so the translations below use it. samplyr’s pps_brewer is Brewer’s general method for any sample size, whereas PPS_BREWER in SAS and SPSS selects two units per stratum.

SAS PROC SURVEYSELECT

proc surveyselect data=frame method=pps_sampford n=50 seed=12345;
  strata region;
  cluster school;
  size enrollment;
run;

sampling_design() |>
  stratify_by(region) |>
  cluster_by(school) |>
  draw(n = 50, method = "pps_sampford", mos = enrollment) |>
  execute(frame, seed = 1)

SAS Allocation with Bounds

proc surveyselect data=frame method=srs n=500 seed=42;
  strata region / alloc=neyman var=variance_data allocmin=2 allocmax=100;
run;

sampling_design() |>
  stratify_by(region, alloc = "neyman", variance = variance_data) |>
  draw(n = 500, min_n = 2, max_n = 100) |>
  execute(frame, seed = 2)

SAS Rounding Control

proc surveyselect data=frame method=sys samprate=0.02 seed=2 round=nearest;
  strata State;
run;

sampling_design() |>
  stratify_by(State) |>
  draw(frac = 0.02, method = "systematic", round = "nearest") |>
  execute(frame, seed = 3)

SPSS CSPLAN

CSPLAN SAMPLE
  /PLAN FILE='myplan.csplan'
  /DESIGN STRATA=region CLUSTER=school
  /METHOD TYPE=PPS_SAMPFORD
  /SIZE VALUE=50
  /MOS VARIABLE=enrollment.

sampling_design() |>
  stratify_by(region) |>
  cluster_by(school) |>
  draw(n = 50, method = "pps_sampford", mos = enrollment) |>
  execute(frame, seed = 4)