Skip to contents

Plan a Sample with svyplan

Source: vignettes/survey-planning.Rmd
survey-planning.Rmd

From a precision target to a draw

Planning answers how many units to issue and how to distribute them. Selection then applies that plan to a frame. svyplan owns the first calculation and samplyr owns the second.

We use the synthetic ken_enterprises frame. Suppose the main estimate is the proportion of exporting enterprises, with a margin of error of three percentage points at 95 percent confidence.

library(samplyr)
library(svyplan)
library(dplyr)

data(ken_enterprises)

p_export <- mean(ken_enterprises$exporter)
target <- n_prop(p = p_export, moe = 0.03)
target
#> Sample size for proportion (wald)
#> n = 504 (p = 0.14, moe = 0.030, deff = 1)
#> expected cases = 68.9

The result is a svyplan_n object. draw() reads its planned sample size directly.

sample <- sampling_design(title = "Enterprise survey") |>
  draw(n = target) |>
  execute(ken_enterprises, seed = 1)

nrow(sample)
#> [1] 504

Use a plausible planning value from earlier rounds or a pilot when the true population value is unknown. For a proportion, p = 0.5 is the conservative choice under the ordinary Wald variance because it maximizes p * (1 - p).

Means, response, and design effects

For a mean, supply the population variance or a prior estimate of it.

revenue_target <- n_mean(
  var = var(ken_enterprises$revenue_millions),
  moe = 30
)
revenue_target
#> Sample size for mean
#> n = 611 (var = 143091.44, moe = 30.000, deff = 1)

Expected nonresponse increases the number issued. Here the target remains the same, but only 85 percent of issued units are expected to respond.

response_target <- n_prop(
  p = p_export,
  moe = 0.03,
  resp_rate = 0.85
)

c(
  full_response = as.integer(target),
  issued_at_85_percent = as.integer(response_target)
)
#>        full_response issued_at_85_percent 
#>                  504                  593

deff expresses the variance relative to simple random sampling under the planning model. It can reflect clustering, unequal weights, or another design feature, but the value should match the estimand being planned.

n_prop(p = p_export, moe = 0.03, deff = 1.5, resp_rate = 0.85)
#> Sample size for proportion (wald)
#> n = 890 gross (net: 756) (p = 0.14, moe = 0.030, deff = 1.50, resp_rate = 0.85)
#> expected cases = 103.3

A cross-sectional design effect does not automatically describe a change or a pooled longitudinal estimate. Correlation across occasions can help one estimand and hurt another.

Read the calculation in reverse

The prec_*() functions evaluate a sample size rather than solving for it.

prec_prop(p = p_export, n = 600)
#> Sampling precision for proportion (wald)
#> n = 600
#> se = 0.0140, moe = 0.0275, cv = 0.1026, rmoe = 0.2011
#> expected cases = 82.0
prec_mean(
  var = var(ken_enterprises$revenue_millions),
  n = 300,
  mu = mean(ken_enterprises$revenue_millions)
)
#> Sampling precision for mean
#> n = 300
#> se = 21.8397, moe = 42.8050, cv = 0.2485, rmoe = 0.4870

This is useful when budget fixes the take. confint() gives the corresponding planning interval, while predict() can compare supported assumptions such as response rates and design effects.

plan <- n_prop(p = p_export, moe = 0.03, deff = 1.5)

predict(
  plan,
  expand.grid(
    deff = c(1, 1.5, 2),
    resp_rate = c(0.8, 0.9, 1)
  )
)
#>   deff resp_rate         n        se  moe        cv      rmoe
#> 1  1.0       0.8  629.7685 0.0153064 0.03 0.1119441 0.2194065
#> 2  1.5       0.8  944.6527 0.0153064 0.03 0.1119441 0.2194065
#> 3  2.0       0.8 1259.5370 0.0153064 0.03 0.1119441 0.2194065
#> 4  1.0       0.9  559.7942 0.0153064 0.03 0.1119441 0.2194065
#> 5  1.5       0.9  839.6913 0.0153064 0.03 0.1119441 0.2194065
#> 6  2.0       0.9 1119.5884 0.0153064 0.03 0.1119441 0.2194065
#> 7  1.0       1.0  503.8148 0.0153064 0.03 0.1119441 0.2194065
#> 8  1.5       1.0  755.7222 0.0153064 0.03 0.1119441 0.2194065
#> 9  2.0       1.0 1007.6296 0.0153064 0.03 0.1119441 0.2194065

Several indicators or domains

n_multi() finds the binding requirement across several means and proportions.

indicators <- data.frame(
  name = c("exporter_rate", "mean_revenue"),
  p = c(p_export, NA),
  var = c(NA, var(ken_enterprises$revenue_millions)),
  moe = c(0.03, 30)
)

multi_target <- n_multi(indicators)
multi_target
#> Multi-indicator sample size
#> n = 611 (binding: mean_revenue)
#> 
#>  name          .n  .cv_target .cv_achieved .binding
#>  exporter_rate 504 0.1119441  0.1016724            
#>  mean_revenue  611        NA         NA    *

A domain-specific plan needs enough observations within every reporting domain. When n_multi() returns a domain table, the sampling design needs a matching stratify_by() so those counts have an operational meaning.

domain_indicators <- data.frame(
  name = "exporter_rate",
  size_class = c("Small", "Medium", "Large"),
  p = c(0.10, 0.22, 0.45),
  cv = 0.10
)

domain_target <- n_multi(domain_indicators, domains = "size_class")

domain_sample <- sampling_design() |>
  stratify_by(size_class) |>
  draw(n = domain_target) |>
  execute(ken_enterprises, seed = 2)

count(domain_sample, size_class, name = "n_sampled")
#> # A tibble: 3 × 2
#>   size_class n_sampled
#>   <fct>          <int>
#> 1 Small            900
#> 2 Medium           355
#> 3 Large            123

Allocate a total across strata

When strata are already defined, n_alloc() distributes a total using the chosen planning criterion. This example uses Neyman allocation for mean revenue, with at least 20 enterprises in every sector.

allocation_frame <- ken_enterprises |>
  summarise(
    N = n(),
    sd = sd(revenue_millions),
    .by = sector
  ) |>
  transmute(stratum = sector, N, sd)

allocation <- n_alloc(
  allocation_frame,
  n = 600,
  alloc = "neyman",
  min_n_stratum = 20
)

as.data.frame(allocation)[, c("stratum", "N", "n_int", "weight")]
#>                         stratum    N n_int   weight
#> 1 Chemicals & Chemical Products  282    26 10.95380
#> 2                  Construction 2343    77 30.54684
#> 3                          Food  856    62 13.79354
#> 4        Hotels and Restaurants 3065    60 50.63273
#> 5           Other Manufacturing 2113   118 17.91983
#> 6                Other Services 6187   158 39.08919
#> 7                        Retail 2158    99 21.84907

The result maps its stratum names to the values used by stratify_by().

allocated_sample <- sampling_design(title = "Allocated enterprise survey") |>
  stratify_by(sector) |>
  draw(n = allocation) |>
  execute(ken_enterprises, seed = 3)

count(allocated_sample, sector, name = "n_sampled")
#> # A tibble: 7 × 2
#>   sector                        n_sampled
#>   <fct>                             <int>
#> 1 Food                                 62
#> 2 Chemicals & Chemical Products        26
#> 3 Other Manufacturing                 118
#> 4 Construction                         77
#> 5 Retail                               99
#> 6 Hotels and Restaurants               60
#> 7 Other Services                      158

Do not also request an allocation in stratify_by(). Either svyplan supplies the stratum takes or samplyr computes them, not both.

Neyman allocation minimizes the planned variance for the variable represented by sd, subject to the stated constraints. A cost-weighted allocation adds the declared costs. Neither choice is automatically optimal for other variables or domains.

Allocate for several indicators at once

Given per-stratum moments for several indicators and a precision target for each, n_alloc() finds the cheapest allocation that satisfies all of them together. measures carries the moments, one row per stratum and indicator, and targets carries one requirement per indicator.

sector_moments <- ken_enterprises |>
  summarise(
    N = n(),
    mean = mean(revenue_millions),
    sd = sd(revenue_millions),
    p = mean(exporter),
    .by = sector
  )

joint_measures <- bind_rows(
  mutate(sector_moments, stratum = sector, name = "revenue", mean, sd, .keep = "none"),
  mutate(sector_moments, stratum = sector, name = "exporter", p, .keep = "none")
)

joint_plan <- n_alloc(
  mutate(sector_moments, stratum = sector, N, .keep = "none"),
  measures = joint_measures,
  targets = data.frame(
    name = c("revenue", "exporter"),
    cv = c(0.10, 0.05)
  )
)

joint_plan
#> Joint constrained allocation (Bethel)
#> question: cheapest design meeting every precision target
#> field design: n = 1997, cost = 1997 (2 targets, all pass)
#> continuous optimum: n = 1996.895, cost = 1997 (integerizing costs +0.01%)
#> binding: exporter@.overall:cv (target 0.05, achieved 0.05)

One requirement determines the sample size and the others have slack. The result reports which.

joint_plan$constraints[, c("name", ".target", ".achieved", ".binding")]
#>       name .target .achieved .binding
#> 1  revenue    0.10 0.0885843    FALSE
#> 2 exporter    0.05 0.0500000     TRUE

The exporter proportion binds. Mean revenue was asked for a CV of 0.10 and the design delivers about 0.089, because the sample the exporter target requires is larger than the one revenue needs. A total sample size on its own would not show that.

Check a plan against the design that runs it

A plan states the precision it expects. Executing the design repeatedly against a known population shows the precision it delivers. Comparing the two is the most direct check that the variance model, the allocation, and the executed design agree with each other.

joint_design <- sampling_design(title = "Enterprise indicators") |>
  stratify_by(sector) |>
  draw(n = joint_plan)

estimate_once <- function(seed) {
  smp <- execute(joint_design, ken_enterprises, seed = seed)
  c(
    revenue = sum(smp$.weight * smp$revenue_millions) / sum(smp$.weight),
    exporter = sum(smp$.weight * smp$exporter) / sum(smp$.weight)
  )
}

replicates <- 200
estimates <- vapply(seq_len(replicates), estimate_once, numeric(2))

Fix the acceptance band before reading the result. The Monte Carlo error on an estimated CV is about 1 / sqrt(2 * (replicates - 1)) in relative terms, so three of those is a defensible band for a check of this size.

truth <- c(
  revenue = mean(ken_enterprises$revenue_millions),
  exporter = mean(ken_enterprises$exporter)
)
tolerance <- 3 / sqrt(2 * (replicates - 1))

check <- data.frame(
  indicator = joint_plan$constraints$name,
  planned_cv = joint_plan$constraints$.achieved,
  empirical_cv = apply(estimates, 1, sd) / abs(rowMeans(estimates)),
  rel_bias = rowMeans(estimates) / truth[joint_plan$constraints$name] - 1,
  row.names = NULL
)
check$ratio <- check$empirical_cv / check$planned_cv
check$within_tolerance <- abs(check$ratio - 1) < tolerance

check
#>   indicator planned_cv empirical_cv      rel_bias     ratio within_tolerance
#> 1   revenue  0.0885843   0.08778954  0.0020154736 0.9910281             TRUE
#> 2  exporter  0.0500000   0.04829564 -0.0007714097 0.9659129             TRUE

Compare each indicator against the precision the plan plans for it, which is .achieved, and not against the target it was asked for. Revenue is checked against 0.089 rather than 0.10 for the reason the previous section gives.

This validates the planning model, not the survey. The frame is used as its own population, so response, coverage, and measurement play no part, and the moments the plan was given are true by construction. It cannot tell you whether those moments describe the population that will actually be surveyed.

Plan a clustered design

Clustering usually increases variance because units in the same cluster are similar. For a one-stage approximation with equal cluster size, the familiar design effect is determined by the homogeneity measure icc and the take per cluster.

design_effect(icc = 0.05, n_per_psu = 10)
#> Planning design effect: 1.4500

n_cluster() can trade the number of clusters against the number of units per cluster under declared costs and homogeneity. Its result is optimal only under that planning model.

cluster_plan <- n_cluster(
  stage_cost = c(500, 50),
  icc = 0.05,
  budget = 30000
)

cluster_plan
#> Optimal 2-stage allocation
#> field design: n_psu = 25 | n_per_psu = 14 -> total n = 350
#> cv = 0.0687, cost = 30000
#> continuous optimum: n_psu = 25.22699 | n_per_psu = 13.78405 (cv = 0.0687, cost = 30000)
#> design df = 24
as.integer(cluster_plan)
#> [1] 25 14

In a two-stage samplyr design, the cluster plan supplies the stage-1 cluster count and the within-cluster take. See ?n_cluster for the feasibility checks and vignette("introduction") for multistage execution.

Field a certainty-aware allocation

A few large PSUs often dominate a stratum. Under proportional-to-size selection their inclusion probability reaches one, so a plan that pretends they are sampled understates what the field work fixes in advance. Given a register with one row per PSU, n_alloc() separates the certainty part from the sampled remainder and sizes both together. Here the ultimate units are households and the register carries each enumeration area’s household count.

ea_register <- data.frame(
  psu_id = c(sprintf("U%02d", 1:8), sprintf("R%02d", 1:8)),
  stratum = rep(c("urban", "rural"), each = 8),
  N = c(600, rep(200, 7), 450, 220, 165, rep(153, 5))
)

certainty_plan <- n_alloc(
  data.frame(stratum = c("urban", "rural"), N = c(2000, 1600), n_per_psu = 10),
  measures = data.frame(
    stratum = c("urban", "rural"), name = "literacy",
    p = c(0.5, 0.4), icc_psu = 0.05
  ),
  targets = data.frame(name = "literacy", cv = 0.10),
  psu = ea_register
)

certainty_plan$detail[, c("stratum", "n_int", "n_psu_certain", "n_psu_draw", "threshold")]
#>   stratum n_int n_psu_certain n_psu_draw threshold
#> 1   urban    86             1          6  238.1445
#> 2   rural    69             1          5  242.4261

The plan’s $psu table classifies every PSU and states the take the operational design fields in it. A certainty PSU carries its own whole take. The remainder carries n_per_psu.

filter(certainty_plan$psu, certainty)
#>   psu_id stratum   N certainty n_take .certainty_source .threshold .distance
#> 1    U01   urban 600      TRUE     26         threshold   238.1445 1.5194783
#> 2    R01   rural 450      TRUE     19         threshold   242.4261 0.8562357

Passing the plan to draw() fields that classification as stored. Stage 1 must be stratified and clustered, with an exact-pik fixed-size PPS method and the register’s own sizes as mos: every certainty PSU enters with probability one, and exactly n_psu_draw more are drawn per stratum. The same plan at stage 2 applies each PSU’s take. The frame here is the household listing the register describes.

households <- ea_register[rep(seq_len(nrow(ea_register)), ea_register$N), ] |>
  mutate(household = row_number(), .by = psu_id)

certainty_sample <- sampling_design(title = "Certainty-aware household survey") |>
  stratify_by(stratum) |>
  cluster_by(psu_id) |>
  draw(n = certainty_plan, method = "pps_systematic", mos = N) |>
  add_stage() |>
  draw(n = certainty_plan) |>
  execute(households, seed = 4)

count(certainty_sample, stratum, name = "n_sampled")
#> # A tibble: 2 × 2
#>   stratum n_sampled
#>   <chr>       <int>
#> 1 rural          69
#> 2 urban          86

The totals equal the plan’s n_int exactly, because both stages field the plan’s own whole-unit numbers. Before anything is drawn, execute() reconciles the frame against the register on PSU identity, stratum, and size, and refuses a frame the plan was not solved for. It also refuses a plan whose remainder draw would push another PSU to probability one, naming the PSU to flag in the register’s certainty column before refitting. The plan’s classification is never silently altered.

Stratifying within PSUs at stage 2 is allowed when an allocation method states how each take splits. A bare stratify_by() is refused there, because a scalar size at a stratified stage means that size in every stratum.

households$sex <- c("f", "m")[households$household %% 2 + 1]

stratified_takes <- sampling_design() |>
  stratify_by(stratum) |>
  cluster_by(psu_id) |>
  draw(n = certainty_plan, method = "pps_systematic", mos = N) |>
  add_stage() |>
  stratify_by(sex, alloc = "proportional") |>
  draw(n = certainty_plan) |>
  execute(households, seed = 4)

count(stratified_takes, stratum, name = "n_sampled")
#> # A tibble: 2 × 2
#>   stratum n_sampled
#>   <chr>       <int>
#> 1 rural          69
#> 2 urban          86

Designs holding a certainty plan serialize with write_design() like any other and replay identically. The file carries the register, so an edited copy is refused at execution rather than fielded.

Plan a repeated survey

A repeated survey often needs precision for a change. That calculation depends on how much of the responding sample appears at both occasions and how strongly the outcome is associated within those shared units.

design_overlap() computes issued-sample overlap from a rotation schedule. A cohort interviewed on four consecutive occasions produces 75 percent overlap between adjacent steady-state samples.

membership <- design_overlap("4")
membership[1]
#> [1] 0.75

At complete response, this issued overlap can be passed directly to n_change(). Below complete response, convert it to the expected overlap among respondents under a stated response-persistence assumption. For illustration, assume response is independent between occasions at rate 0.80. A shared issued unit then responds again with probability 0.80, so respondent overlap is issued overlap multiplied by that rate.

change_response_rate <- 0.80
respondent_overlap <- membership[1] * change_response_rate

change_target <- n_change(
  p = c(p_export, p_export + 0.05),
  moe = 0.02,
  overlap = respondent_overlap,
  overlap_cor = 0.6
)

change_target
#> Sample size for change (proportion scale)
#> n = 1667 per occasion (p = 0.137 to 0.187, moe = 0.02, deff = 1)
#> overlap = 0.6, overlap_cor = 0.6 (64.3% of the independent variance)

overlap_cor is an outcome assumption, not another name for membership overlap. Use prior rounds or a sensitivity range rather than deriving it from the rotation pattern.

Recruit a rotating panel

n_panel() translates a cross-sectional responding target into recruitment under initial response and conditional retention.

panel_plan <- n_panel(
  n_prop(p = p_export, moe = 0.03),
  retention = c(0.90, 0.95, 0.95),
  resp_rate = 0.80,
  design = "rotating",
  start = "immediate"
)

panel_plan
#> Panel recruitment (rotating, 4-wave life)
#> entrants: 177 per occasion -> 505 responding, pooled over 4 cohorts
#> proportion (wald): se = 0.01531, moe = 0.03, cv = 0.112
#> 
#>  wave retention n_resp se      moe    
#>  1              142    0.02891 0.05666
#>  2    0.9       127    0.03047 0.05973
#>  3    0.95      121    0.03126 0.06128
#>  4    0.95      115    0.03208 0.06287
#> 
#> # summary() for the launch, the loss and per-wave cv

A rotating svyplan_panel contains several counts because start-up and later intakes are different operations.

intake_n <- ceiling(panel_plan$n_entrants)

fieldwork_counts <- c(
  startup = panel_plan$n_cohorts * intake_n,
  intake_each_occasion = intake_n,
  expected_respondents = panel_plan$n_resp
)

fieldwork_counts
#>              startup intake_each_occasion expected_respondents 
#>             708.0000             177.0000             503.8148

Use startup for the immediate start-up master and partition it into four cohort ages. Use intake_each_occasion for every later refreshment draw. draw(n = panel_plan) is intentionally unsupported because it cannot know which operation the caller means.

n_panel() treats retention as an unbroken chain. It does not yet turn a gapped life such as 4-8-4 into a recruitment and reserve-cohort plan. It also does not estimate attrition, adjust nonresponse weights, or combine cohort weights. Those are separate design and analysis decisions.

The complete execution handoff is in vignette("rotating-panels").

What planning does not decide

Planning results depend on their inputs. Before treating a number as a field requirement, record:

  • the estimand and reporting domains
  • the source of variances, proportions, and homogeneity assumptions
  • whether the count means issued units or expected respondents
  • the response and retention model
  • the design effect and the estimand it describes
  • the rounding rule and any finite-frame limit

svyplan does not choose a sampling algorithm or validate a register. samplyr does not decide the precision target. Keeping that boundary explicit makes it possible to revisit assumptions without silently changing the selection design.

Next steps

References