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"), which includes a three-stage
household sample. The sections on joins, export and domains apply to any
supported sample.
We use a small stratified EA sample so the analysis examples run
independently. Sample inspection is covered in the “Inspecting the
result” section of vignette("introduction").
sample <- sampling_design(title = "Burkina Faso EA survey") |>
stratify_by(region, alloc = "proportional") |>
draw(n = 300) |>
execute(bfa_eas, seed = 1)variance_estimators() reports, from the design and the
frame and before anything is drawn, which variance estimators the design
supports, which only approximately, and which the export would refuse,
with the condition each would raise. A stratified first stage rules out
survey’s JK1, for instance.
estimators <- variance_estimators(get_design(sample), bfa_eas)
estimators[c("estimator", "status", "class")]
#> # A tibble: 11 × 3
#> estimator status class
#> <chr> <chr> <chr>
#> 1 Taylor linearization supported NA
#> 2 Linearization, joint probabilities not applicable NA
#> 3 Jackknife (JK1) refused samplyr_error_svrep_conve…
#> 4 Stratified jackknife (JKn) supported NA
#> 5 Balanced repeated replication refused samplyr_error_svrep_conve…
#> 6 Fay's balanced repeated replication refused samplyr_error_svrep_conve…
#> 7 Bootstrap supported NA
#> 8 Rescaled bootstrap supported NA
#> 9 Multistage rescaled bootstrap supported NA
#> 10 Rao-Wu-Yue-Beaumont bootstrap supported NA
#> 11 Random groups supported NAJoin 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.
Nonresponse and calibration
Keep every issued row and record each unit’s disposition next to its
outcomes. Missing outcomes do not change the selection weights, and
neither dropping nonrespondents before export nor
na.rm = TRUE is a nonresponse adjustment. The adjustment
and any calibration belong to the analysis. Applied to the exported
design, they leave .weight holding the base weights while
the adjusted weights live in the design object.
set.seed(3)
issued <- sample_with_data |>
mutate(disposition = if_else(runif(n()) < 0.85, "respondent",
"nonrespondent"))
issued |> as.data.frame() |> count(disposition)
#> disposition n
#> 1 nonrespondent 37
#> 2 respondent 263
set.seed(4)
issued_replicates <- as_svrepdesign(issued, type = "rwyb", replicates = 500)
adjusted <- svrep::redistribute_weights(
issued_replicates,
reduce_if = disposition == "nonrespondent",
increase_if = disposition == "respondent",
by = "region"
)
respondents <- subset(adjusted, disposition == "respondent")
#> Warning in `[.svyrep.design`(x, r, ): Fewer observations than degrees of
#> freedom in subsetThe adjustment moves each nonrespondent’s weight to the respondents
of its region, in every replicate as well as in the full-sample weights,
so the replicate variance carries it. It assumes that response is
unrelated to the outcome within a region, which the analyst has to
justify. survey warns that the respondent subset has fewer
rows than the design’s degrees of freedom. The replicate variance is
computed as for any domain and is not affected. Calibration then matches
the respondents to known population totals, here the number of urban EAs
in the frame.
ea_totals <- c(`(Intercept)` = nrow(bfa_eas),
urban_ruralUrban = sum(bfa_eas$urban_rural == "Urban"))
calibrated <- calibrate(respondents, ~urban_rural, population = ea_totals)
svymean(~consumption_per_capita, calibrated)
#> mean SE
#> consumption_per_capita 79988 3726.5weights(calibrated) returns the adjusted and calibrated
weights, while issued$.weight keeps the base weights, so
the two stay apart for reporting. The example uses RWYB replicates
because svrep’s weight adjustments act on a replicate design.
?variance_estimation gives what RWYB assumes for each
selection method.
How each selection method’s variance is estimated is in
?variance_estimation, and ?as_svydesign and
?as_svrepdesign give the 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, pps = FALSE)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 multi-stage 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 without replacement 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.
Systematic sampling uses an SRS-style variance approximation, and
systematic PPS uses Brewer’s. Neither sees frame order, so the direction
of the error depends on it. as_svydesign() warns unless the
caller explicitly accepts or refuses the approximation through
systematic_variance.
Strata with one sampled unit
A stratum that ends up with one sampled unit outside certainty,
because its allocation was one or because certainty units took the rest,
has no within-stratum variance estimate. execute() says so
when it selects such a stratum, and as_svydesign() warns
again at export (samplyr_warning_lonely_psu) while
survey is left at its default
survey.lonely.psu = "fail", under which estimation would
stop. Setting the option, as below, also silences that export
warning.
one_per_region <- sampling_design() |>
stratify_by(region) |>
draw(n = 1) |>
execute(bfa_eas, seed = 5)
#> Stage 1: 13 strata take a single unit outside certainty.
#> ℹ Strata with one unit: "Boucle du Mouhoun", "Cascades", "Centre",
#> "Centre-Est", "Centre-Nord", and 8 more.
#> ℹ A stratum needs two selections for its variance to be estimated. survey stops
#> on a stratum with one unless its strata are collapsed, and a replicate method
#> has nothing to resample in it.
#> ℹ To estimate within each stratum, allocate at least two units outside
#> certainty to it. `min_n` bounds a stratum's whole take, certainty units
#> included.
old <- options(survey.lonely.psu = "adjust")
svytotal(~households, as_svydesign(one_per_region))
#> total SE
#> households 3684462 721634
options(old)"adjust" centers a lonely unit at the population mean,
which is conservative. Collapsing each such stratum with a similar one
before export is the other usual choice, and it is a design decision the
analyst has to make. Allocating at least two units outside certainty to
every stratum avoids the problem at the source. RWYB replicates refuse a
singleton whose variance is needed, except under Poisson sampling and,
with lonely.psu = "certainty", at the final stage.
?variance_estimation gives the certainty case.
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.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 7074466This transformation does not repair a poor variance estimator, and it applies only when the estimate is positive. It helps where the sampling distribution of a positive total is skewed, as for a small domain, and matters little for a full-population total. For a skewed ratio, neither a log-ratio nor a Fieller interval is uniformly reliable, and the ratio estimator itself can carry material finite-sample bias.
Replicate weights and PPS variance
For a supported single-phase design, as_svrepdesign()
creates replicate weights through survey::as.svrepdesign(),
through the optional svrep package when
type = "rwyb" is requested, or from the replicates of
execute(reps = R) when type = "random_groups"
is requested.
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.
?variance_estimation gives how each selection method’s
variance is estimated and how exact it is.
For independent Poisson selection, including clustered and
multi-stage designs, the generic replicate types, "auto"
included, refuse the sample, so 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.6RWYB retains the variation in Poisson sample size. It also supports
SRS, independent draws with replacement, and Brewer, CPS, Sampford and
systematic PPS stages, alone or mixed across 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. Multi-stage 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.
Other designs
The same export functions handle designs this article does not work through, and each help page carries a runnable example.
| Design | Where it is documented |
|---|---|
| A second phase selected from the first-phase sample |
?execute for the selection,
?as_svydesign for the survey::twophase()
export |
| Exact or approximate joint inclusion probabilities for PPS |
?joint_expectation, passed to
as_svydesign(pps = )
|
| Independent samples from frames that overlap on one population |
?stack_frames,
?as_svydesign.frame_stack and
?as_svrepdesign.frame_stack
|
| Weights shared to a linked population | ?share_weights |
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 and the contract separating design weights from estimation
weights.