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 TRUEThese 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
}
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_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 × 9
#> 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.950The 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.
Interpreting the results
For this short validation, practical acceptance bands are:
-
abs(rel_bias_pct) <= 2, -
0.85 <= se_ratio <= 1.15, -
0.90 <= coverage_95 <= 0.98.
The last two are Monte Carlo widths rather than statistical
thresholds, and both work out to about three Monte Carlo standard errors
at R = 199.
At nominal 95% coverage the Monte Carlo standard error of the
coverage estimate is \(\sqrt{0.95 \times 0.05
/ R}\), which is about 1.5 percentage points at
R = 199. The band 0.90 to 0.98 is roughly that either side
of 0.95, a little wider below than above.
The empirical SD in se_ratio is itself estimated from
R draws, with a relative standard error of about \(1 / \sqrt{2(R - 1)}\), which is 5% at
R = 199. The band 0.85 to 1.15 is three of those.
Reading the bands this way is what makes them a judgement rather than
two magic numbers: a se_ratio of 1.12 is inside the band
because a correct design will land there often at this replication
count, not because a 12% error would be acceptable in a real
estimate.
acceptance <- mc_results |>
mutate(
pass_bias = abs(rel_bias_pct) <= 2,
pass_se = se_ratio >= 0.85 & se_ratio <= 1.15,
pass_cov = coverage_95 >= 0.90 & coverage_95 <= 0.98,
pass_all = pass_bias & pass_se & pass_cov
) |>
select(design, pass_bias, pass_se, pass_cov, pass_all)
acceptance
#> # A tibble: 3 × 5
#> design pass_bias pass_se pass_cov pass_all
#> <chr> <lgl> <lgl> <lgl> <lgl>
#> 1 SRSWOR TRUE TRUE TRUE TRUE
#> 2 Stratified SRS TRUE TRUE TRUE TRUE
#> 3 Two-stage PPS->SRS TRUE TRUE TRUE TRUE
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.0154A pass across these three designs says that:
- the deterministic invariants on
samplyrdesign metadata hold exactly, which is a proof for the sample drawn and not a Monte Carlo result - nothing in the exported
surveydesigns is far enough from correct design-based inference forR = 199to detect.
The second is the weaker claim, and it is the one that matters here. The acceptance bands are about three Monte Carlo standard errors wide, so a bias or a variance error smaller than that will pass unnoticed. A pass is consistent with correctness rather than evidence for it.
R is fixed at 199 as a practical default so the vignette
builds quickly. Rerunning with a larger R narrows every
band above in proportion to \(1/\sqrt{R}\) and is what to do when the
question is how large a residual error the run can rule out.