From selection to estimation
The usual workflow has four parts. Select the sample, collect
outcomes, join them to the selected units, and export the result to
survey or srvyr.
samplyr records the strata, clusters, weights, and
finite population corrections during selection. A safe one-to-one join
preserves that metadata. An operation that adds, removes, or duplicates
sampled rows changes the realization and prevents export.
Inspect the sample before fieldwork
We begin with a stratified sample of enumeration areas.
data(bfa_eas)
sample <- sampling_design(title = "Burkina Faso EA survey") |>
stratify_by(region, alloc = "proportional") |>
draw(n = 300) |>
execute(bfa_eas, seed = 1)
summary(sample)
#> ── Sample Summary: Burkina Faso EA survey ──────────────────────────────────────
#>
#> ℹ n = 300 of 44,570 | stages = 1/1 | seed = 1
#>
#> ── 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 300frame_summary() exposes the same realization at pool
level.
frame_summary(sample, detail = "pool") |>
select(region, N, n_target, n_expected, n_realized, take_rate)
#> # A tibble: 13 × 6
#> region N n_target n_expected n_realized take_rate
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Boucle du Mouhoun 5009 34 34 34 0.00679
#> 2 Cascades 2508 17 17 17 0.00678
#> 3 Centre 3888 26 26 26 0.00669
#> 4 Centre-Est 2941 20 20 20 0.00680
#> 5 Centre-Nord 3402 23 23 23 0.00676
#> 6 Centre-Ouest 3723 25 25 25 0.00672
#> 7 Centre-Sud 1612 11 11 11 0.00682
#> 8 Est 5505 37 37 37 0.00672
#> 9 Hauts-Bassins 4839 32 32 32 0.00661
#> 10 Nord 2930 20 20 20 0.00683
#> 11 Plateau-Central 1662 11 11 11 0.00662
#> 12 Sahel 4144 28 28 28 0.00676
#> 13 Sud-Ouest 2407 16 16 16 0.00665For this feasible fixed-size design, the target, expectation, and
realized take coincide. Under Bernoulli or Poisson selection the
realized count is random. Under unequal-probability selection,
take_rate describes the realized pool count and is not a
unit-level inclusion probability.
Join collected outcomes
The frame need not contain the survey outcomes. Join one collected record to each sampled unit using a stable identifier.
library(dplyr)
set.seed(2)
collected <- tibble(
ea_id = sample$ea_id,
consumption_per_capita = round(rlnorm(
nrow(sample),
meanlog = 11,
sdlog = 0.65
)),
below_threshold = rbinom(nrow(sample), size = 1, prob = 0.25)
)
sample_with_data <- sample |>
left_join(collected, by = "ea_id")
class(sample_with_data)
#> [1] "tbl_sample" "tbl_df" "tbl" "data.frame"The result remains a tbl_sample because the join is one
to one and leaves the selected rows intact. Check key uniqueness before
joining when collected data can contain repeated visits, roster members,
or duplicate submissions.
Export to survey
as_svydesign() maps the recorded design to a
survey object.
library(survey)
svy <- as_svydesign(sample_with_data)
svy
#> Stratified Independent Sampling design
#> svydesign(ids = ~1, strata = ~region, weights = ~.weight, fpc = ~.fpc_1,
#> data = data, nest = TRUE)The collected variables can now be analyzed with standard
survey functions.
svymean(~consumption_per_capita, svy)
#> mean SE
#> consumption_per_capita 78450 3411.2
svymean(~below_threshold, svy)
#> mean SE
#> below_threshold 0.26342 0.0256
svytotal(~population, svy)
#> total SE
#> population 20024627 788500For a multistage sample, every represented stage contributes its sampling-unit and finite-population terms. This preserves the stage structure. It does not turn an approximate stage variance treatment into an exact one.
Fixed-size PPS methods use Brewer’s approximation by default. In that
phrase, Brewer names the survey variance approximation and
not necessarily the selection algorithm. The approximation’s accuracy
depends on the population and design. Exact joint inclusion
probabilities can be supplied for methods that provide them.
Equal-probability systematic sampling uses an SRS-style variance
approximation. Its direction depends on frame ordering.
as_svydesign() warns unless the caller explicitly accepts
or refuses that approximation through
systematic_variance.
Estimate a domain correctly
Convert the complete sample first, then subset the survey design.
urban_svy <- subset(svy, urban_rural == "Urban")
svymean(~consumption_per_capita, urban_svy)
#> mean SE
#> consumption_per_capita 80326 7961.3Filtering the tbl_sample first treats the observed
domain sample as fixed. The point estimate can agree, but the variance
estimate is generally wrong and is often too small because domain
membership is random under the design (Särndal, Swensson, and Wretman 1992,
ch. 5.8).
urban_rows <- sample_with_data |>
filter(urban_rural == "Urban")
as_svydesign(urban_rows)
#> Error in `as_svydesign()`:
#> ! `as_svydesign()` requires a sample that still matches its executed
#> design.
#> ✖ Rows were removed, added, or duplicated after `execute()`.
#> ℹ For domain (subpopulation) analysis, convert the full sample first, then
#> subset the design: `subset(as_svydesign(full_sample), condition)`, or with
#> srvyr: `as_survey_design(full_sample) |> filter(condition)`.
#> ℹ To subsample an executed sample, run a second phase: `sampling_design() |>
#> draw(...) |> execute(full_sample)`.One-to-one joins, row reordering, and new ordinary columns remain
supported. Changing the row set or an internal column such as
.weight or .fpc_1 does not.
With srvyr, use the same order. Convert first and then
filter the survey object.
library(srvyr)
#>
#> Attaching package: 'srvyr'
#> The following object is masked from 'package:stats':
#>
#> filter
sample_with_data |>
as_survey_design() |>
filter(urban_rural == "Urban") |>
summarise(
mean_consumption = survey_mean(consumption_per_capita)
)
#> # A tibble: 1 × 2
#> mean_consumption mean_consumption_se
#> <dbl> <dbl>
#> 1 80326. 7961.Replicate weights and PPS variance
For a supported single-phase design, as_svrepdesign()
creates replicate weights through
survey::as.svrepdesign().
replicate_design <- as_svrepdesign(sample, type = "auto")
svymean(~households, replicate_design)
#> mean SE
#> households 65.744 2.7501Replicate methods are not automatic proofs of exactness. Generic PPS
bootstraps do not recreate spatial selection, hard cube constraints, or
a random Poisson sample size. as_svydesign() and
as_svrepdesign() report or refuse combinations that cannot
be represented under their documented contracts. See their reference
pages before selecting a variance route for an advanced design.
For supported PPS methods, joint_expectation()
reconstructs exact or approximate second-order quantities according to
the method.
sampford_sample <- sampling_design() |>
stratify_by(region) |>
cluster_by(ea_id) |>
draw(n = 5, method = "pps_sampford", mos = households) |>
execute(bfa_eas, seed = 2025)
joints <- joint_expectation(sampford_sample, bfa_eas, stages = 1)
exact_pps_design <- as_svydesign(
sampford_sample,
pps = ppsmat(joints[[1]])
)
svymean(~households, exact_pps_design)
#> mean SE
#> households 76.339 7.6495The Sampford matrix above is exact. CPS, systematic PPS, and Poisson
selection also have exact second-order calculations in the supported
cases. Generalized Brewer, SPS, Pareto, and unconstrained cube use the
documented high-entropy approximation. Spatial methods and bounded cube
do not provide a joint matrix. ?joint_expectation is the
authoritative method table.
A user-supplied pps matrix describes one stage. For a
multistage design, supplying it exports only the represented PPS stage.
Omit it when the intended analysis is the full multistage linearization
with the default stage-level variance treatment.
Frames that overlap on one population
Two registers of the same population, each incomplete: a landline frame and a cell frame, an area frame and a list frame. The samples are selected independently, one design and one seed each, and everything the feature contains happens after that.
stack_frames() records which frames each sampled unit
belongs to. Every component has to carry every frame’s membership
column, not only its own, because a unit selected from frame A has to
say whether it was also listed in frame B, and that is the information
the compositing needs.
population <- data.frame(
person_id = 1:600,
spend = round(stats::rlnorm(600, log(120), 0.4), 2),
in_landline = rep(c(TRUE, FALSE), times = c(400, 200)),
in_cell = rep(c(FALSE, TRUE), times = c(150, 450))
)
frames <- stack_frames(
landline = sampling_design() |>
draw(n = 80) |>
execute(population[population$in_landline, ], seed = 1),
cell = sampling_design() |>
draw(n = 120) |>
execute(population[population$in_cell, ], seed = 2),
membership = c(landline = "in_landline", cell = "in_cell"),
key = person_id
)
summary(frames)
#> ── Frame Stack Summary ─────────────────────────────────────────────────────────
#>
#> ℹ 2 frames | 200 rows | key person_id
#>
#> Frames
#> • landline: 80 rows | in_landline | seed 1
#> • cell: 120 rows | in_cell | seed 2
#>
#> Domains
#> • cell: 49 rows
#> • cell+landline: 117 rows
#> • landline: 34 rows
#>
#> ℹ A unit listed in two frames is counted once per frame that selected it.The domain counts are rows, not distinct people. Someone listed in both frames and selected from both appears twice, once per frame that selected them, and that is exactly the double count a composite weight has to resolve.
Two export routes
as_svydesign() hands the stack to
survey::multiframe(). Each component is exported on its
own, so it keeps its strata, clusters and variance treatment, and the
compositing sits on top.
svytotal(~spend, as_svydesign(frames))
#> total SE
#> spend 80964 2958.4theta = NULL, the default, is the multiplicity
estimator: a unit reached by both frames contributes half its weight
through each. An explicit theta is Hartley’s constant
factor (Hartley
1962), and it belongs to the first frame of
the stack, so reversing the components and wanting the same estimator
means asking for 1 - theta.
svytotal(~spend, as_svydesign(frames, theta = 0.74))
#> total SE
#> spend 79841 3065.5samplyr never forwards a missing theta to
survey. Left unset there, it defaults to the ratio of the
frames’ mean sampling weights, which is a data-dependent stand-in for
frame size over sample size rather than a neutral choice.
The replicate route builds one system per frame and combines them in blocks: in a replicate column belonging to one frame, only that frame varies and every other stays at its full-sample weight. The combined variance is then the sum of the frames’ own contributions, which is what independent selection from each frame gives.
svytotal(~spend, as_svrepdesign(frames, type = "bootstrap", replicates = 100))
#> total SE
#> spend 80964 2852.7The two routes are not interchangeable, and the differences are the reason to know which one you are on.
as_svydesign() |
as_svrepdesign() |
|
|---|---|---|
| Number of frames | Two. survey::multiframe() composites two |
Any number |
| Compositing | Multiplicity, or Hartley with a stated theta
|
The same, with theta above two frames refused |
| A component with shared weights | Refused | Supported |
| Replicate methods | Not applicable | May differ between frames; type = "auto" picks one per
component |
The two-frame limit belongs to the export and is stated there, not at
the design layer: stack_frames() records any number,
exactly as execute() chains any number of phases while the
two-phase export stops at two.
Where the theory comes from
The multiplicity estimator for overlapping frames is the generalized
weight share method with the standardized link matrix, because the
number of links to a unit is the number of frames containing it (Lohr 2021, sec. 3.2; Mecatti 2007). That is why a stack
composed of frames and a sample whose weights were shared through a link
table are handled by one mechanism rather than two, and why
vignette("design-semantics") documents them together.
Analysis checklist
- Inspect the selected counts and weights before fieldwork.
- Verify that collected-data keys are unique at the sampled-unit level.
- Join outcomes without changing the selected row set.
- Convert the complete sample to a survey design.
- Define domains on the converted design.
- State every variance approximation used by an advanced selection method.
- Where the estimates combine frames, say which compositing factor was
used and which frame an explicit
thetawas applied to.
Use vignette("serialization") to archive the selection
and vignette("design-semantics") for the detailed export
contracts, the weight share method, and the contract separating design
weights from estimation weights.