Analyze a samplyr Sample

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.

The main path reaches this article after vignette("introduction") and vignette("three-stage-sampling"). The sections on joins, export and domains apply to an ordinary sample. Read the later two-phase or overlapping-frame examples when your design uses them.

We use a small stratified EA sample so the analysis examples run independently. Sample inspection is covered in the selection tutorials.

sample <- sampling_design(title = "Burkina Faso EA survey") |>
  stratify_by(region, alloc = "proportional") |>
  draw(n = 300) |>
  execute(bfa_eas, seed = 1)

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.

Keep the issued rows when some units do not respond: use a left join and record a response indicator. Missing outcomes do not change the selection weights. Neither dropping those rows before export nor using na.rm = TRUE is, by itself, a nonresponse adjustment. Choosing response adjustments or calibration belongs to the downstream analysis and its assumptions.

The compact selection and inference table is in ?selection-methods. ?as_svydesign and ?as_svrepdesign give the detailed export restrictions.

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 788500

For 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.3

Filtering 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.

For a strictly positive domain total, a log-scale delta interval is an available finite-sample sensitivity analysis. It keeps the endpoints positive and allows asymmetric uncertainty.

urban_total <- svytotal(~population, urban_svy)
urban_estimate <- unname(coef(urban_total)[[1]])
urban_se <- unname(SE(urban_total)[[1]])
z <- qnorm(0.975)

exp(log(urban_estimate) + c(-1, 1) * z * urban_se / urban_estimate)
#> [1] 4780414 7074466

This transformation does not repair a poor variance estimator. It applies only when the estimate is positive. Extended PPS qualification found a clear coverage improvement for positive domain totals, but not for full-population totals. Log-ratio and Fieller intervals did not provide a uniformly reliable correction for skewed ratios, whose point estimators can also have material finite-sample bias.

Replicate weights and PPS variance

For a supported single-phase design, as_svrepdesign() creates replicate weights through survey::as.svrepdesign(), or through the optional svrep package when type = "rwyb" is requested.

replicate_design <- as_svrepdesign(sample, type = "auto")
svymean(~households, replicate_design)
#>              mean     SE
#> households 65.744 2.7501

Replicate 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 independent Poisson selection, including clustered and multistage designs, request Rao-Wu-Yue-Beaumont replicates explicitly:

poisson_sample <- sampling_design() |>
  draw(frac = 0.4, method = "bernoulli") |>
  execute(data.frame(id = 1:100, y = 1:100), seed = 12)
set.seed(42)
poisson_replicates <- as_svrepdesign(
  poisson_sample, type = "rwyb", replicates = 1000
)
svytotal(~y, poisson_replicates)
#>    total    SE
#> y 4597.5 683.6

RWYB retains the variation in Poisson sample size that generic fixed-size bootstraps can lose. It also supports mixed SRS, independent draws with replacement, and fixed-size PPS WOR stages. PPS WOR variance is an approximation and is reported as such. Increase the replicate count to reduce simulation error. svrep’s estimate_boot_sim_cv() helps assess it for chosen estimates. No new sampling verbs are needed.

Two-phase replicate export is not yet supported. Multistage RWYB export requires every selected parent to remain represented in the final sample. Later Poisson stages require a complete frame digest to check this. A noncertainty singleton is refused when its variance is not estimable, while a singleton Poisson sample can contribute a valid variance estimate.

When exporting a materialized wave or a second-phase sample with as_svydesign(), current analysis columns take precedence over same-named columns in the retained first-phase sample. Missing current measurements remain missing. Keep earlier measurements under separate names if you need both. Identifiers and sampling metadata retain their design meanings.

Two phases: screening followed by a detailed interview

A second phase samples units already selected at phase 1, using information collected in that phase. It is a new selection law on the realized sample. A later stage in a household design instead samples children within selected parents. Continuing that design with a household listing is not a new phase.

The bundled ken_enterprises is entirely synthetic. Here a cheap screening interview determines phase-two strata, and only the subsample receives a detailed interview. The stable enterprise_id links the phases. At phase 1, cluster_by(enterprise_id) declares that identity. Each such unit has one row.

screening_design <- sampling_design("Establishment screening") |>
  cluster_by(enterprise_id) |>
  draw(n = 500)

phase1 <- execute(screening_design, ken_enterprises, seed = 71)
set.seed(72)
phase1 <- phase1 |>
  mutate(
    revenue_baseline = revenue_millions,
    screen_score = log1p(employees) + rnorm(n(), sd = 0.5),
    priority = if_else(screen_score >= 4, "high", "ordinary")
  )
phase1 |> as.data.frame() |> count(priority)
#>   priority   n
#> 1     high  83
#> 2 ordinary 417

The synthetic screening variables are observed for every phase-one unit, including those that will not enter phase 2. They remain in the retained phase-one sample so the two-phase export can reconstruct the selection pools.

interview_design <- sampling_design("Detailed establishment interview") |>
  stratify_by(priority) |>
  draw(n = c(high = 90, ordinary = 60))

phase2 <- execute(interview_design, phase1, seed = 73)
#> Warning: Stage 1: sample size exceeded the pool population in 1 of 2 pools.
#> ✖ Requested 150 units, selected 143.
#> ℹ Capped pools: "high".
#> ℹ Inspect with `frame_summary(sample, detail = "pool")` and the capped column.

set.seed(74)
phase2 <- phase2 |>
  mutate(
    revenue_millions = revenue_baseline * exp(rnorm(n(), 0, 0.15)),
    audit_hours = rgamma(n(), shape = 3, rate = 0.5)
  )

audit_hours is a new phase-two measurement. revenue_millions is an updated measurement with the same name as its baseline counterpart. Export uses the current phase-two values. revenue_baseline keeps the earlier values available explicitly. All detailed interviews are complete in this example.

two_phase_svy <- as_svydesign(phase2)
svytotal(~revenue_millions + audit_hours, two_phase_svy)
#>                    total     SE
#> revenue_millions 1212828 233654
#> audit_hours       103231   8010

# The export uses the current observations and the combined phase weights.
stopifnot(isTRUE(all.equal(
  unname(coef(svytotal(~revenue_millions, two_phase_svy))),
  sum(phase2$.weight * phase2$revenue_millions)
)))

These estimates illustrate inference for the synthetic establishment population, not substantive estimates for Kenya. The example uses SRS at phase 1 and stratified SRS at phase 2, a supported two-phase path. PPS, WR and Poisson phase designs have additional export restrictions. Phase chains longer than two and two-phase replicate export are currently unsupported. See ?as_svydesign before changing the selection methods.

Joint information for PPS variance

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.6495

The 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.

An exact systematic PPS matrix can contain zero pair probabilities. Exact design quantities do not then guarantee an estimable or stable sample variance through survey::ppsmat(). A sampled matrix cannot show whether unseen population pairs have zero probability, so as_svydesign() warns when a ppsmat object is supplied for systematic PPS. It does not refuse the route, because full pair positivity can hold at high sampling fractions.

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 76409 2925.1

theta = 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 74818 2917.5

The expected multiframe estimator instead needs each unit’s selection probability in every frame containing it. Supply the original registers through overlaps = exante_overlaps(list(area = area_register, list = list_register), by = c(person_id = "person_id")) when building the stack, using the stack’s component names and key, then export with estimator = "expected". The resolver requires exact inclusion probabilities by default. SPS and Pareto targets require allow_approximate = TRUE in exante_overlaps(). This accepts an approximation that can affect design unbiasedness. The stack records quality by frame, and exported designs retain it in their samplyr_overlap_probability_quality attribute. For inspecting one register, exante_probabilities() returns the same quality beside each probability.

samplyr 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 76409 3036.1

The 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

1. Inspect the selected counts and weights before fieldwork.
2. Verify that collected-data keys are unique at the sampled-unit level.
3. Join outcomes without changing the selected row set.
4. Convert the complete sample to a survey design.
5. Define domains on the converted design.
6. State every variance approximation used by an advanced selection method.
7. Where the estimates combine frames, say which compositing factor was used and which frame an explicit theta was 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.

References

Hartley, H. O. 1962. “Multiple Frame Surveys.” In Proceedings of the Social Statistics Section, American Statistical Association, 203–6.

Lohr, Sharon L. 2021. “Multiple-Frame Surveys for a Multiple-Data-Source World.” Survey Methodology 47 (2): 229–63.

Mecatti, Fulvia. 2007. “A Single Frame Multiplicity Estimator for Multiple Frame Surveys.” Survey Methodology 33 (2): 151–57.

Särndal, Carl-Erik, Bengt Swensson, and Jan Wretman. 1992. Model Assisted Survey Sampling. Springer.