Validation on Synthetic Populations • samplyr

Why this vignette

This vignette runs a set of smoke tests for samplyr on synthetic finite populations with known truths. It uses two complementary layers:

Deterministic checks on design metadata and invariants.
Monte Carlo checks for estimator bias, standard-error calibration, and 95% CI coverage.

The goal is not to benchmark speed, and it is not a proof of correctness. A run at R = 199 can detect a design that is plainly wrong. It cannot establish that one is right. A pass says that the package behaved consistently with correctness on these three designs at this replication count. That is the level of assurance a smoke test provides. The scale of the Monte Carlo uncertainty is reported below so the bands can be read for what they are.

Synthetic population with known properties

library(dplyr)
set.seed(20260302)

make_population <- function(
    n_strata = 6L,
    psu_per_stratum = 50L,
    units_per_psu = 30L) {
  n_psu <- n_strata * psu_per_stratum
  stratum <- rep(LETTERS[seq_len(n_strata)], each = psu_per_stratum * units_per_psu)
  psu <- rep(sprintf("PSU%03d", seq_len(n_psu)), each = units_per_psu)

  psu_mos <- round(runif(n_psu, min = 80, max = 420))
  psu_re <- rnorm(n_psu, sd = 18)
  stratum_re <- rep(
    seq(-12, 16, length.out = n_strata),
    each = psu_per_stratum * units_per_psu
  )

  x <- rep(psu_mos, each = units_per_psu) + rnorm(n_psu * units_per_psu, sd = 10)
  y <- 200 + stratum_re + rep(psu_re, each = units_per_psu) + 0.08 * x + rnorm(n_psu * units_per_psu, sd = 20)

  tibble(
    id = seq_len(n_psu * units_per_psu),
    stratum = factor(stratum),
    psu = factor(psu),
    psu_mos = rep(psu_mos, each = units_per_psu),
    x = x,
    y = y
  )
}

population <- make_population()

truth <- population |>
  summarise(
    N = dplyr::n(),
    Y_total = sum(y),
    Y_mean = mean(y)
  )

truth
#> # A tibble: 1 × 3
#>       N  Y_total Y_mean
#>   <int>    <dbl>  <dbl>
#> 1  9000 2008269.   223.

Layer 1: deterministic checks

library(samplyr)
library(survey)

srs_sample <- sampling_design() |>
  draw(n = 400) |>
  execute(population, seed = 1001)

strat_sample <- sampling_design() |>
  stratify_by(stratum, alloc = "proportional") |>
  draw(n = 400) |>
  execute(population, seed = 1002)

pps_two_stage <- sampling_design() |>
  add_stage("PSU") |>
    stratify_by(stratum) |>
    cluster_by(psu) |>
    draw(
      n = 8,
      method = "pps_brewer",
      mos = psu_mos,
      certainty_size = 390
    ) |>
  add_stage("SSU") |>
    draw(n = 6) |>
  execute(population, seed = 1003)

stratum_sizes <- population |>
  count(stratum, name = "N_h")

strat_weight_by_stratum <- strat_sample |>
  summarise(weight_sum = sum(.weight), .by = stratum) |>
  left_join(stratum_sizes, by = "stratum")

# all(logical(0)) is TRUE, so the certainty check has to establish that the
# set it quantifies over is non-empty before it means anything.
n_certainty <- sum(pps_two_stage$.certainty_1)

checks <- tibble(
  check = c(
    "SRS: sum(weights) = N",
    "Stratified SRS: sum_h weights_h = N_h",
    "Two-stage PPS: .weight = .weight_1 * .weight_2",
    "Two-stage PPS: the certainty_size threshold selected some PSUs",
    "Two-stage PPS certainty: certainty PSUs have .weight_1 = 1"
  ),
  pass = c(
    isTRUE(all.equal(sum(srs_sample$.weight), truth$N)),
    isTRUE(all(abs(strat_weight_by_stratum$weight_sum - strat_weight_by_stratum$N_h) < 1e-8)),
    isTRUE(all.equal(pps_two_stage$.weight, pps_two_stage$.weight_1 * pps_two_stage$.weight_2, tolerance = 1e-10)),
    n_certainty > 0,
    n_certainty > 0 && all(pps_two_stage$.weight_1[pps_two_stage$.certainty_1] == 1)
  )
)

checks
#> # A tibble: 5 × 2
#>   check                                                          pass 
#>   <chr>                                                          <lgl>
#> 1 SRS: sum(weights) = N                                          TRUE 
#> 2 Stratified SRS: sum_h weights_h = N_h                          TRUE 
#> 3 Two-stage PPS: .weight = .weight_1 * .weight_2                 TRUE 
#> 4 Two-stage PPS: the certainty_size threshold selected some PSUs TRUE 
#> 5 Two-stage PPS certainty: certainty PSUs have .weight_1 = 1     TRUE

These are structural, non-random checks. If they fail, the design object is wrong before analysis begins.

The fourth row is not decoration. all() over an empty set returns TRUE, so a certainty check written without it would pass on a design in which certainty_size = 390 happened to select nothing at all, and the pass would say nothing about certainty handling.

Layer 2: Monte Carlo validation

Below we repeatedly sample from the same finite population and analyze with survey. For each design we report:

rel_bias_pct: relative bias of \(\hat{\bar{Y}}\),
bias_z: the same bias divided by its own Monte Carlo standard error,
se_ratio: empirical SD of estimates divided by mean estimated SE,
coverage_95: fraction of 95% CIs covering the true finite-population mean.

rel_bias_pct and bias_z measure the same thing on different scales, and only the second is comparable across designs. The true mean here is near 200, so a 2% relative bias is about 4 units, which is large next to the Monte Carlo noise: relative bias against a large mean is a generous criterion. bias_z divides by the Monte Carlo standard error of the estimated bias, so a value near zero says the bias is not distinguishable from simulation noise, and a value past 2 or 3 says it is, whatever the units of y happen to be.

This vignette fixes Monte Carlo replication at R = 199.

evaluate_design <- function(design_fn, population, true_mean, R = 199L, seed_init = 3000L) {
  est <- numeric(R)
  se <- numeric(R)
  cover <- logical(R)

  for (i in seq_len(R)) {
    sample_i <- design_fn(population, seed_init + i)
    svy_i <- as_svydesign(sample_i)
    fit <- svymean(~y, svy_i)

    est[i] <- unname(coef(fit)[1])
    se[i] <- unname(SE(fit)[1])

    ci <- suppressWarnings(confint(fit, level = 0.95))
    cover[i] <- ci[1, 1] <= true_mean && ci[1, 2] >= true_mean
  }

  # Delete one independent simulation at a time. This accounts for both
  # parts of the SD/mean-SE ratio and their dependence within each draw.
  leave_one_ratio <- vapply(seq_len(R), function(i) {
    sd(est[-i]) / mean(se[-i])
  }, numeric(1))
  ratio_mcse <- sqrt((R - 1) * mean(
    (leave_one_ratio - mean(leave_one_ratio))^2
  ))
  coverage_ci <- binom.test(sum(cover), R)$conf.int

  tibble(R = R,
         bias = mean(est - true_mean),
         rel_bias_pct = 100 * mean(est - true_mean) / true_mean,
         bias_z = mean(est - true_mean) / (sd(est) / sqrt(R)),
         emp_sd = sd(est),
         mean_se = mean(se),
         se_ratio = sd(est) / mean(se),
         coverage_95 = mean(cover),
         mcse_rel_bias_pct = 100 * sd(est) / sqrt(R) / abs(true_mean),
         mcse_se_ratio = ratio_mcse,
         coverage_lower = coverage_ci[1],
         coverage_upper = coverage_ci[2])
}

mcse_coverage <- function(R, p = 0.95) {
  sqrt(p * (1 - p) / R)
}

design_srs <- function(pop, seed) {
  sampling_design() |>
    draw(n = 400, method = "srswor") |>
    execute(pop, seed = seed)
}

design_strat <- function(pop, seed) {
  sampling_design() |>
    stratify_by(stratum, alloc = "proportional") |>
    draw(n = 400, method = "srswor") |>
    execute(pop, seed = seed)
}

design_two_stage <- function(pop, seed) {
  sampling_design() |>
    add_stage("PSU") |>
      stratify_by(stratum) |>
      cluster_by(psu) |>
      draw(n = 8, method = "pps_brewer", mos = psu_mos) |>
    add_stage("SSU") |>
      draw(n = 6, method = "srswor") |>
    execute(pop, seed = seed)
}

R_mc <- 199L
R_mc
#> [1] 199

mc_results <- bind_rows(
  evaluate_design(design_srs, population, truth$Y_mean, R = R_mc) |>
    mutate(design = "SRSWOR"),
  evaluate_design(design_strat, population, truth$Y_mean, R = R_mc) |>
    mutate(design = "Stratified SRS"),
  evaluate_design(design_two_stage, population, truth$Y_mean, R = R_mc) |>
    mutate(design = "Two-stage PPS->SRS")
) |>
  select(design, everything())

mc_results
#> # A tibble: 3 × 13
#>   design       R    bias rel_bias_pct bias_z emp_sd mean_se se_ratio coverage_95
#>   <chr>    <int>   <dbl>        <dbl>  <dbl>  <dbl>   <dbl>    <dbl>       <dbl>
#> 1 SRSWOR     199 -0.0344      -0.0154 -0.344   1.41    1.41    0.998       0.950
#> 2 Stratif…   199 -0.128       -0.0573 -1.34    1.34    1.32    1.01        0.940
#> 3 Two-sta…   199  0.0816       0.0366  0.355   3.24    2.99    1.09        0.950
#> # ℹ 4 more variables: mcse_rel_bias_pct <dbl>, mcse_se_ratio <dbl>,
#> #   coverage_lower <dbl>, coverage_upper <dbl>

The third design’s first stage is pps_brewer, and as_svydesign() exports it with Brewer’s approximation to the joint inclusion probabilities. Its mean_se is therefore an approximate quantity, and se_ratio for that row compares the empirical spread against an approximation rather than against an exact variance. The first two rows have no such caveat.

Separate tolerances from Monte Carlo uncertainty

A practical tolerance states how much error is acceptable for an application. Monte Carlo uncertainty states how precisely this simulation measured the error. They answer different questions. For illustration, choose tolerances of two percentage points of relative bias, 15% around an SD/mean-SE ratio of one, and three percentage points around nominal 95% coverage. These are teaching choices, not universal standards for qualifying a survey design.

The table below compares those fixed tolerances with 95% Monte Carlo intervals. Bias uses the empirical SD divided by the square root of R. The ratio uses a jackknife over independent simulation runs, accounting for uncertainty in both its numerator and denominator. These two intervals are normal approximations. Coverage uses an exact binomial interval. The intervals measure simulation uncertainty, not uncertainty in a particular survey estimate.

precision <- bind_rows(
  mc_results |>
    transmute(design, metric = "Relative bias (%)", estimate = rel_bias_pct,
              lower = estimate - 1.96 * mcse_rel_bias_pct,
              upper = estimate + 1.96 * mcse_rel_bias_pct,
              acceptable_lower = -2, acceptable_upper = 2),
  mc_results |>
    transmute(design, metric = "SD / mean SE", estimate = se_ratio,
              lower = pmax(0, estimate - 1.96 * mcse_se_ratio),
              upper = estimate + 1.96 * mcse_se_ratio,
              acceptable_lower = 0.85, acceptable_upper = 1.15),
  mc_results |>
    transmute(design, metric = "95% coverage", estimate = coverage_95,
              lower = coverage_lower, upper = coverage_upper,
              acceptable_lower = 0.92, acceptable_upper = 0.98)
) |>
  mutate(assessment = case_when(
    lower >= acceptable_lower & upper <= acceptable_upper ~ "Within tolerance",
    upper < acceptable_lower | lower > acceptable_upper ~ "Outside tolerance",
    TRUE ~ "Inconclusive at this replication count"
  ))

precision |>
  select(design, metric, estimate, lower, upper, assessment) |>
  knitr::kable(digits = 3,
               col.names = c("Design", "Metric", "Estimate", "MC lower",
                             "MC upper", "Assessment"))

Design	Metric	Estimate	MC lower	MC upper	Assessment
SRSWOR	Relative bias (%)	-0.015	-0.103	0.072	Within tolerance
Stratified SRS	Relative bias (%)	-0.057	-0.141	0.026	Within tolerance
Two-stage PPS->SRS	Relative bias (%)	0.037	-0.165	0.239	Within tolerance
SRSWOR	SD / mean SE	0.998	0.903	1.092	Within tolerance
Stratified SRS	SD / mean SE	1.014	0.918	1.110	Within tolerance
Two-stage PPS->SRS	SD / mean SE	1.085	1.004	1.167	Inconclusive at this replication count
SRSWOR	95% coverage	0.950	0.910	0.976	Inconclusive at this replication count
Stratified SRS	95% coverage	0.940	0.897	0.968	Inconclusive at this replication count
Two-stage PPS->SRS	95% coverage	0.950	0.910	0.976	Inconclusive at this replication count

An estimate can fall inside a tolerance while its Monte Carlo interval extends outside. That result is inconclusive, not a pass. Conversely, a small bias can be statistically distinguishable from zero (bias_z) while remaining practically tolerable. Report both size and uncertainty.

tibble(
  R = R_mc,
  mcse_coverage_95 = mcse_coverage(R_mc, p = 0.95)
)
#> # A tibble: 1 × 2
#>       R mcse_coverage_95
#>   <int>            <dbl>
#> 1   199           0.0154

At R = 199, nominal 95% coverage has Monte Carlo SE of about 1.5 percentage points. Resolving a narrow coverage tolerance may therefore need substantially more draws. Increasing R typically narrows the Monte Carlo intervals at roughly the inverse-square-root rate. It does not shrink the chosen practical tolerances. The interval calculations above are recomputed from the new simulations. The tolerance constants remain unchanged unless the application’s requirements change.

These are pointwise diagnostics for three particular populations/designs, not simultaneous certification of the package. Deterministic checks establish the stated identities for the realized examples. The simulations assess mean bias, the SD/mean-SE ratio and interval coverage in these cases. They do not establish unbiased variance estimation for all totals or all populations. Random-size sampling, extreme MOS, small pools and other compositions need separate checks. See ?selection-methods for the inference overview and ?as_svydesign and ?as_svrepdesign for the detailed export restrictions.