This vignette documents what each samplyr operation
assumes, what it computes, and what is preserved or lost at
survey export. The contracts stated here are asserted.
vignette("validation") checks them against synthetic
populations with known truths and reports what a run at that replication
count can and cannot establish.
Assumptions at a glance
The following table summarizes what samplyr assumes and
what the user must verify. Each row is expanded in the sections
below.
| Assumption | Applies to | samplyr behavior | User must verify |
|---|---|---|---|
| Strata are exhaustive and disjoint | stratify_by() |
Partitions the frame by unique combinations of strata variables | Strata variables define a valid partition |
| Independent sampling across strata | stratify_by() |> draw() |
Each stratum sampled with a separate call to
sondage
|
Correct by construction when strata are disjoint |
| Cluster variables constant within clusters | cluster_by() |
Validated at execution: MOS and strata variables may not vary within a cluster | Cluster ID uniquely identifies the group |
| Conditional independence across stages | add_stage() |
Stage \(k\) sampling depends only on the set of units selected at stage \(k - 1\) | Standard multi-stage assumption that fails if a later frame depends on earlier randomization |
| FPC stored as population count | execute() |
.fpc_k contains \(N_h\) or a cluster count and is converted
to \(\pi_i\) at survey export for PPS
stages |
Frame size at execution is treated as the population, one of several definitions when a population changes between waves |
| Weights are \(1/\pi_i\), unadjusted | All methods | No calibration, trimming, or non-response adjustment | Downstream adjustment is the analyst’s responsibility |
| PRN values are iid \(\text{U}(0, 1)\) | draw(prn = ...) |
Validated: strictly in \((0, 1)\), no NAs, numeric | Independence and uniformity of the PRN column |
Notation
Throughout this vignette:
- \(U = \{1, \ldots, N\}\): finite population (the frame).
- \(S \subseteq U\): selected sample.
- \(\pi_i = \Pr(i \in S)\): first-order inclusion probability.
- \(\pi_{kl} = \Pr(k \in S \text{ and } l \in S)\): joint (second-order) inclusion probability.
- \(w_i = 1 / \pi_i\): design weight (sampling weight before any adjustment).
- \(N_h\): population size of stratum \(h\) (\(h = 1, \ldots, H\)).
- \(n_h\): sample size drawn from stratum \(h\).
- \(K\): number of sampling stages.
- \(\pi_i^{(k \mid S^{(k-1)})}\): conditional inclusion probability of unit \(i\) at stage \(k\), given the set \(S^{(k-1)}\) selected at stage \(k - 1\).
execute() produces design weights only.
Calibration, trimming, and non-response adjustment are out of scope and
belong in the estimation layer (survey, srvyr,
or similar).
One declared transformation departs from that, and says so in the
object it returns: share_weights() produces
estimation weights for a population other than the one
that was sampled. Its semantics, and the contract that keeps the two
kinds of weight apart, are the subject of a section below.
What each verb means
sampling_design(): frame-independent plan
A sampling_design object is a specification, not a
computation. Column names referenced in stratify_by(),
cluster_by(), and draw() are stored as strings
(deferred resolution). No data is touched until execute()
binds the design to a frame. This separation means the same design can
be applied to different frames or re-executed with different seeds.
design <- sampling_design(title = "SRS WOR") |>
stratify_by(region) |>
draw(n = 50, method = "srswor")
# No data has been sampled yet
design
#> ── Sampling Design: SRS WOR ────────────────────────────────────────────────────
#>
#> ℹ 1 stage
#>
#> ── Stage 1 ─────────────────────────────────────────────────────────────────────
#> • Strata: region
#> • Draw: n = 50 (per stratum), method = srswor
stratify_by(): partition into independent
sub-populations
stratify_by() partitions the frame into \(H\) strata defined by the Cartesian product
of its variables:
\[ U = U_1 \cup U_2 \cup \cdots \cup U_H, \qquad U_h \cap U_{h'} = \emptyset \; \text{for } h \neq h'. \]
Sampling within each stratum is independent. This holds by
construction because samplyr invokes sondage
algorithms separately for each stratum.
When an allocation method is specified, the total sample size \(n\) is distributed across strata:
| Allocation | Formula | Requires |
|---|---|---|
"equal" |
\(n_h = n / H\) | |
"proportional" |
\(n_h \propto N_h\) | |
"neyman" |
\(n_h \propto N_h \, S_h\) | variance |
"optimal" |
\(n_h \propto N_h \, S_h / \sqrt{C_h}\) |
variance, cost
|
"power" |
\(n_h \propto \text{CV}_h \cdot X_h^q\) |
cv, importance
|
Each formula gives a factor \(w_h\). The allocation is the solution of
\[ n_h = \operatorname{clamp}(\lambda w_h, \; \ell_h, \; u_h), \qquad \ell_h = \min(\texttt{min\_n}, N_h), \qquad u_h = \min(\texttt{max\_n}, N_h) \]
with \(\lambda\) chosen so that
\(\sum_h n_h = n\). The population size
\(N_h\) is an upper bound whether or
not max_n was supplied, so a fixed-size without-replacement
design satisfies \(0 \le n_h \le N_h\)
and \(\sum_h n_h = n\) for every
feasible request.
Saturated strata are fixed at their bound and removed from the problem, then \(\lambda\) is recomputed over the strata still free. Redistribution therefore follows the same factors \(w_h\) the method is defined by, which is what keeps a Neyman allocation Neyman after a stratum saturates: redistributing in proportion to unused capacity instead would silently change the criterion. When every free stratum has \(w_h = 0\) the criterion cannot distinguish them and the remainder is split equally in stratum order.
One consequence is worth stating plainly: a named rule can be
departed from when a stratum is too small to hold its share. With
populations \((10, 490, 500)\) and
\(n = 300\), "equal"
yields 10/145/145, not 100/100/100. execute() reports this
once per run.
The constrained real-valued targets are then rounded with largest-remainder (Hare–Niemeyer) rounding that preserves \(\sum n_h = n\) and stays inside every bound.
With-replacement and Poisson-multinomial designs draw units
repeatedly, so \(N_h\) is not a bound
on the number of draws and only min_n and
max_n apply. Per-stratum sizes supplied directly, rather
than through alloc, are instructions and are never
redistributed.
cluster_by(): define the sampling unit
cluster_by() declares the variable(s) that identify the
primary sampling unit (PSU) at the current stage. Operationally,
samplyr:
- Groups the frame by the cluster variable(s).
- Draws a sample of clusters (not individual units).
- Returns all population units within selected clusters.
Any variable used in PPS size measures (mos) or in
stratification must be constant within each cluster.
samplyr validates this at execution time and raises an
error if cluster-level variables vary within groups.
cluster_by() is purely structural. It does not affect
inclusion probabilities directly. They are determined by
draw() applied to the cluster-level aggregates.
draw(): specify the selection mechanism
draw() attaches a selection method and sample size (or
fraction) to the current stage. The inclusion probabilities and weights
depend on the method class:
Equal-probability without replacement
(srswor, systematic,
bernoulli):
\[ \pi_i = \frac{n}{N} \quad (\text{or } \pi_i = f \text{ for bernoulli}), \qquad w_i = \frac{1}{\pi_i}. \]
Fixed-size for srswor and systematic and
random-size for bernoulli. Systematic sampling induces
implicit stratification by row order. Use control to
specify the sort variable.
PPS methods (pps_brewer,
pps_systematic, pps_cps,
pps_sampford, pps_poisson,
pps_sps, pps_pareto):
\[ \pi_i = n \, \frac{\text{mos}_i}{\sum_j \text{mos}_j} \]
computed by sondage::inclusion_prob(), which iteratively
caps units with \(\pi_i \geq 1\) and
redistributes. These “certainty” units (\(\pi_i = 1\)) receive stage weight 1 and
contribute zero first-stage variance. In multi-stage designs, final
.weight still compounds all stages, so certainty at one
stage does not imply overall weight 1. For Poisson sampling, \(\pi_i = f \cdot N \cdot \text{mos}_i / \sum
\text{mos}\), capped at 1 (random-size design).
Two of the seven attain that \(\pi_i\) only approximately.
pps_sps and pps_pareto are order sampling
designs, which fix the sample size at the cost of the marginals: the
quantity above is the target inclusion probability
handed to the algorithm, and the realized inclusion probability is close
to it but not equal, approaching it as \(n\) and \(N\) grow (Tillé 2020, sec. 5.11; Rosén 1997). Both sit in the
"approximate" probabilities tier. Everywhere below that
reads \(\pi_i\) off
.weight inherits that distinction. The other five methods
satisfy it exactly.
With-replacement and PMR (srswr,
pps_multinomial, pps_chromy):
The quantity of interest is the expected number of hits \(E(n_i) = n \, \text{mos}_i / \sum
\text{mos}\). Each draw is an independent selection. If unit
\(i\) is selected \(k\) times, samplyr creates
\(k\) rows with a .draw
column indexing each hit. The weight per draw is \(w_i = 1 / p_i\), where \(p_i\) is the single-draw selection
probability. No finite population correction applies (FPC = \(\infty\)).
Chromy’s sequential PPS (pps_chromy) is classified as
probability minimum replacement (PMR): each unit receives \(\lfloor E(n_i) \rfloor\) or \(\lceil E(n_i) \rceil\) hits. When all \(E(n_i) < 1\), this reduces to WOR.
Balanced and spatially balanced sampling
(cube, lpm2, scps):
The cube method (Deville and Tillé 2004) selects
a fixed-size sample that satisfies (or nearly satisfies) the balancing
equations \(\sum_{i \in S} a_{ij} / \pi_i =
\sum_{i \in U} a_{ij}\) for auxiliary variables \(a_j\) specified via aux.
Without mos, \(\pi_i =
n/N\) and probabilities are equal. With mos, \(\pi_i\) is computed by
sondage::inclusion_prob(). A bound() marker
adds hard adjacent-integer category-count constraints.
When stratified, the stratified cube algorithm (Chauvet 2009) runs flight phases per stratum, then a global flight phase, then per-stratum landing phases.
LPM2 and SCPS instead use the coordinates supplied through
spread. They do not accept cube auxiliaries or count
bounds. At most 2 stages may use a balanced-family method. PRN and
certainty selection are not supported.
PRN-compatible methods (bernoulli,
pps_poisson, pps_sps,
pps_pareto):
When prn is specified, the permanent random numbers are
passed through to sondage. Selection becomes deterministic
given the PRN values, so the seed argument to
execute() has no effect on the selection outcome.
add_stage(): delimit stages
add_stage() is syntactic: it closes the current stage
and opens a new one. It carries no statistical content of its own.
Stages are numbered \(1, 2, \ldots, K\)
in definition order.
execute(): materialize the sample
execute() binds the design to one or more frames and
runs the selection:
- For each stage \(k\), the frame is restricted to units belonging to clusters selected at stage \(k - 1\) (or the full frame for \(k = 1\)).
- Strata and cluster groupings are resolved against the (subsetted) frame.
- The selection algorithm from
draw()is invoked viasondage. - The
.weight,.fpc,.draw, and.certaintycolumns are computed and suffixed with the stage number.
One frame or one per stage. The number of frames
schedules the stages. execute(design, hierarchy) runs every
stage against one table, and
execute(design, schools, classes, students) gives each
stage its own register, mapped by position. Any other count is refused.
The frames may equally be held as a value with
execute(design, list(schools, classes, students)). A name
on a list element is a diagnostic label for that frame.
A register is supplied whole. samplyr restricts it to the units its parent stage selected and carries onto it the variables earlier stages introduced, so a class register needs neither pre-filtering to the sampled schools nor a copy of the school-level stratification column. A register that omits such a carried variable is fine. One whose own copy disagrees with the parent is an error.
Partial execution.
execute(design, psu_frame, stages = 1) runs only the first
stage and returns a tbl_sample. That object stores both the
unchanged design and the realized PSU selection. Use it as the first
argument of the later call and supply the new listing as the frame:
execute(psu_sample, listing_frame). Execution resumes at
the next unexecuted stage, carries forward the PSU weights and FPC, and
remains one multi-stage design.
The stages argument may select any initial contiguous
batch beginning at stage 1 and on a partial tbl_sample, it
may select the next contiguous batch. Omitting it executes all remaining
stages, with one exception. If more than one stage remains and a single
frame is supplied, that frame could be the next stage’s register or a
hierarchy covering the rest. The two draw different samples, so
execute() asks for stages rather than
guessing.
The three spellings are one implementation. A
three-stage design executed against one hierarchy, against three
registers in a single call, or stage by stage as a continuation runs the
same stage transition each time. All three draw the same sample,
provided the continuation runs under a single RNG stream and names
stages on every intermediate call:
one_call <- execute(design, schools, classes, students, seed = 7)
hierarchical <- execute(design, hierarchy, seed = 7)
continued <- withr::with_seed(7, {
s1 <- execute(design, schools, stages = 1)
s2 <- execute(s1, classes, stages = 2)
execute(s2, students, stages = 3)
})Note where the seed goes. One seed wraps the whole chain, and the individual calls take none. Seeding each call separately is a different and equally valid RNG boundary, but it is not this equivalence and will select different units. The equivalence is a tested invariant rather than an incidental property. The continued sample is therefore the one the single call would have drawn.
Two-phase sampling. A new phase starts with a new
design: execute(design2, phase1_sample). Here the phase-1
tbl_sample is the frame argument, not the execution object.
samplyr records a previous-phase link and compounds the phase-1 and
phase-2 weights multiplicatively. The conditional phase-2 weight is
computed at survey export time. Thus
execute(psu_sample, listing_frame) means “continue this
design’s stages,” whereas execute(new_design, prior_sample)
means “sample a new phase.”
Weight compounding across stages
Single-stage. The design weight is:
\[ w_i = \frac{1}{\pi_i}. \]
Multi-stage. Let \(\pi_i^{(k)}\) denote the conditional inclusion probability of unit \(i\) at stage \(k\), given the set of clusters selected at all prior stages. The compound weight is:
\[ w_i = \prod_{k=1}^{K} \frac{1}{\pi_i^{(k)}} = \prod_{k=1}^{K} w_i^{(k)}. \]
samplyr computes this by joining each stage’s weight to
the previous stage’s weight on the shared cluster variables and
multiplying. The per-stage weights are preserved as
.weight_1, .weight_2, etc. The
.weight column always equals their product.
two_stage <- sampling_design() |>
add_stage("Districts") |>
stratify_by(province) |>
cluster_by(district) |>
draw(n = 2, method = "srswor") |>
add_stage("EAs") |>
draw(n = 5, method = "srswor") |>
execute(zwe_eas, seed = 123456)
#> Warning: Stage 1: sample size exceeded the pool population in 1 of 10 pools.
#> ✖ Requested 20 units, selected 19.
#> ℹ Capped pools: "Bulawayo".
#> ℹ Inspect with `frame_summary(sample, detail = "pool")` and the capped column.
# Compound weight equals the product of per-stage weights
all.equal(two_stage$.weight,
two_stage$.weight_1 * two_stage$.weight_2)
#> [1] TRUEThe warning is the design meeting the frame. Bulawayo is a
single-district province, so a take of two districts selects the one
that exists. execute() reports every pool it could not
fill, once per stage, because a capped pool makes the design
non-self-weighting. The weights above still compound exactly and
Bulawayo’s first-stage weight is simply 1.
Two-phase. When a new phase-2 design is executed
with the phase-1 tbl_sample as its frame, the phase-1
weight is stored internally and multiplied into the final
.weight. The conditional phase-2 weight \(w_{\text{cond}} = w_{\text{overall}} /
w_{\text{phase 1}}\) is computed during survey export for use
with survey::twophase().
Finite population corrections
samplyr stores the FPC as a population count in
.fpc_k. The interpretation depends on the method, and the
conversion for survey export happens automatically in
as_svydesign():
| Method class |
.fpc_k stores |
Passed to survey as | survey interpretation |
|---|---|---|---|
| EP-WOR (srswor, systematic) | \(N_h\) (stratum pop.) | \(N_h\) directly | Sampling fraction \(f_h = n_h / N_h\) |
| PPS-WOR (pps_brewer, …) | \(N_h\) (stratum pop.) | \(1 / w_k = \pi_i\) | Inclusion probability |
| WR / PMR (srswr, chromy, …) | \(\infty\) (synthetic) | \(\infty\) | No FPC (Hansen–Hurwitz variance) |
The PPS conversion merits explanation.
survey::svydesign() interprets FPC values less than 1 as
inclusion probabilities and values \(\geq
1\) as population sizes. samplyr stores \(N_h\) uniformly for all methods, then at
export time, as_svydesign() creates a per-row column
.fpc_pi_k = 1 / .weight_k = \pi_i for PPS-WOR stages. The
transformation is lossless: the per-stage weight already encodes \(\pi_i\), so nothing is approximated in
converting it. What is approximate for pps_sps and
pps_pareto is the quantity being transformed, since their
\(\pi_i\) is a target rather than the
realized inclusion probability.
For WR and PMR methods, a synthetic column filled with
Inf is created. The survey package interprets this as “no
finite population correction,” which gives the Hansen–Hurwitz variance
estimator.
The population the correction refers to.
samplyr treats the frame as executed as the population, so
\(N_h\) is counted from the frame the
stage drew from. That is one convention among several. Where the
sampling fraction is not negligible and the population itself changes
between occasions, the finite population correction matters more and
admits several competing definitions, particularly for measures of
change across coordinated or overlapping samples (Tam 1984; Laniel 1988; Nordberg 2000; Berger 2004b).
samplyr does not choose among them for you: it records the
frame size it saw, and a different convention has to be imposed
downstream.
Independence across stages
samplyr assumes that the stage-\(k\) selection mechanism depends on stage
\(k - 1\) only through the
set of selected clusters, not through the specific
randomization that produced it. Formally:
\[ \Pr(S^{(k)} \mid S^{(k-1)}) \text{ depends on } S^{(k-1)} \text{ as a set.} \]
This is the standard assumption in multi-stage sampling theory (Cochran 1977, ch. 11; Särndal, Swensson, and Wretman 1992, ch. 4.3).
Consequence for variance estimation. Under this
assumption, the recursive multi-stage variance formula applies (Särndal, Swensson, and
Wretman 1992, ch. 4.3). For methods with a supported
linearization estimator, as_svydesign() exports every
executed stage to survey::svydesign(): one ids
term, one fpc term, and (when the stage is stratified) one
strata term per stage. A first-stage census, for example,
correctly attributes all variance to the later stages.
A stage that samples with replacement is handled differently: later stages nested in it contribute nothing beyond the between-draw variability, which the Hansen–Hurwitz estimator already captures.
Bounded cube, LPM2, and SCPS do not currently have supported linearization estimators. They use the generic bootstrap approximation described below.
Weights for a population other than the one sampled
Two situations produce a weight for units that were never in the frame. Both are transformations of an executed sample and neither changes the draw.
Links to a second population
share_weights() is the generalized weight share method
(Lavallée
2007). A sample is selected from population \(U^A\), and estimates are wanted for a
second population \(U^B\) linked to it:
children through their parents, establishments through their
enterprises, persons through their dwellings. Every target unit \(k\) in a reached cluster \(i\) receives
\[ w_i = \sum_j \frac{I(j \in S^A)}{\pi_j} \frac{L_{ji}}{L_i}, \]
the source weights carried across the links and divided by the population number of links to that cluster. Assigning one weight per cluster is what makes unit-level and cluster-level estimates of the same total agree.
Two properties of that formula decide the whole API.
The first is that \(\pi_j\) is
needed only for units actually selected (Lavallée 2007, sec.
2.2.2), which is exactly what .weight already
holds. That is why this is a post-execute transformation rather than a
sixth verb.
The second is that \(L_i\) sums over
the whole of \(U^A\), not over
the sample. Counting a supplied link table therefore asserts
that the table is a complete population register, and getting that wrong
understates every weight whose links are only partly recorded. The
assertion has no default: either name a column holding the population
multiplicity, or state it with complete_links().
dwellings <- data.frame(dwelling_id = 1:200)
dwelling_sample <- sampling_design() |>
draw(n = 20) |>
execute(dwellings, seed = 1)
people <- data.frame(
person_id = 1:400,
dwelling_id = rep(1:200, each = 2)
)
links <- data.frame(
dwelling_id = rep(1:200, each = 2),
person_id = 1:400
)
people_sample <- share_weights(
dwelling_sample,
targets = people,
links = links,
by = c(dwelling_id = "dwelling_id"),
to = c(person_id = "person_id"),
within = dwelling_id,
multiplicity = complete_links()
)
# Each person is linked to the one dwelling they live in, so the target
# weight is the source weight and the totals agree.
c(
dwellings = sum(dwelling_sample$.weight),
people = sum(people_sample$.weight)
)
#> dwellings people
#> 200 400The target cluster is never inferred, because the three ways of
setting it produce three different weights: a bare column is the
ordinary clustered method, NULL makes every target unit its
own cluster, and extend_links() eliminates clusters by
extending the links across them (Lavallée 2007, sec. 5.3).
Links need not be counted. weighted_links() replaces the
0/1 indicator with a non-negative importance, which costs no theory as
long as each target cluster totals more than zero (Lavallée 2007, sec.
4.5). Counting is the case where every link counts for one,
and the two are one implementation rather than two.
The weight contract
A shared weight is an estimation weight. It is not \(1/\pi\) for anything, because the target units have no inclusion probabilities of their own, and the design the object still carries describes the selection of the source population. So the result records which kind of weight it holds, and every statistical consumer states what it does about that record rather than treating the sample as ordinary.
| Consumer | On a shared weight |
|---|---|
design_effect(), effective_n()
|
Supported. They read .weight and nothing else |
as_svrepdesign() |
Supported. It replicates the source design |
| The data-manipulation methods | Preserve the record |
as_svydesign() |
Supported. It exports the source-target contributions |
joint_expectation(), varcomp()
|
Refuse |
A further execute(), and the panel verbs |
Refuse |
Each refusal has its own condition class and all of them are also
catchable as samplyr_error_weight_contract. Neither
as_tbl_sample() nor an ordinary dplyr reconstruction can
launder an estimation weight back into a design-weight sample.
Where the coverage finding fires
A target cluster with no link to the source population can never be
surveyed, and the estimator understates totals by exactly its share.
share_weights() records that finding and does
not warn about it, which departs on purpose from the
convention that a condition fires where it arises.
The reason is the composed case. A transformation cannot know whether
its result will stand alone or become one component of a design over
several frames, and there a cluster one frame cannot reach is the entire
reason the second frame exists. So the warning fires where the estimate
is formed: as_svrepdesign() reports it for a single
transformation, and for a stack_frames() collection it is
evaluated over the union of the frames and fires once, naming only the
clusters no frame reaches.
One component’s silence about a cluster counts as coverage only if that component was describing the same target clusters. Where that cannot be established the union is reported as not established, rather than guessed.
There are two variance routes, and each is correct for a different reason.
as_svrepdesign() replicates the source
design and applies the recorded operator inside every replicate, so the
target rows are never resampled. Sharing after replication is a
different, and wrong, number.
as_svydesign() exports the source-target
contributions: one row per link, weighted by the recorded
coefficient times the source unit’s design weight. The identity behind
it is that the weight share total is the Horvitz-Thompson total of a
variable derived on the source units,
\[ \sum_i w_i y_i = \sum_j \frac{I(j \in S)}{\pi_j} z_j, \qquad z_j = \sum_i \frac{L_{ji}}{L_i} y_i, \]
and z depends on whatever is being analyzed, so it
cannot be formed in advance. Expanding the contributions lets
survey form it inside each sampling unit instead. That is
exact for any link structure, with no condition on how many source units
reach a target.
An unequal-probability or random-size source design is refused there. Both take their variance from a structure indexed by the rows of the source sample, and those rows stop being the sampled units once each appears once per contribution.
Frames that overlap on one population
The second situation is two registers of the same population, each
incomplete. stack_frames() collects samples selected
independently from them, and the theory turns out to be the same one:
the multiplicity estimator for overlapping frames is
the weight share method with the standardized link matrix, since the
number of links to a unit is the number of frames containing it (Lohr 2021, sec. 3.2; Mecatti 2007).
What a stack records is which frames each sampled unit belongs to.
What it does not do is composite the weights, because
which compositing factor is appropriate depends on the estimand and on
the design effects, so it is an estimation-time choice.
theta = NULL at export is the multiplicity estimator,
giving every frame reaching a unit an equal share; an explicit
theta is Hartley’s constant factor (Hartley 1962)
and belongs to the first frame of the stack.
A stack also does not assert that the frames together cover the
target population. That is an assumption about the registers, not
something the call can establish, which is why the verb is not named
after their union. See vignette("survey-analysis") for the
two export routes and what each one can carry.
What is lost at survey export
as_svydesign() translates a tbl_sample into
a survey::svydesign() object. The translation is faithful
for standard variance estimators, but some samplyr metadata
is not representable in the survey framework.
1. Element sampling before later stages
An unclustered element-sampling stage followed by further stages is
not nested cluster sampling: the later selections operate on the
realized element sample, not within larger units, so there is no
hierarchy for the multi-stage variance recursion (or for
replicate-weight resampling) to use. This is phase sampling.
as_svydesign() raises an error for this shape and the
supported path is a two-phase sample: execute the element stage under
its first-phase design, execute a new second-phase design with that
sample as its frame, and export with as_svydesign(), which
uses survey::twophase(). A final unclustered stage
is fully supported: it is exported with a synthesized row-identity
ids term and its own FPC.
2. Joint inclusion probabilities
By default, as_svydesign() uses Brewer’s variance
approximation (Berger
2004a) for fixed-size PPS WOR methods, including Sampford.
Here Brewer names the variance estimator, not the selection algorithm.
The approximation derives \(\pi_{kl}\)
from marginal \(\pi_i\) without
computing the full \(N \times N\)
matrix.
joint_expectation() computes second-order quantities
directly instead. What it returns depends on the method:
- Exact for CPS, Sampford, systematic, and Poisson.
- The high-entropy approximation (Hájek 1964; Brewer and Donadio 2003), which is \(O(N^2)\), for generalized Brewer, SPS, Pareto, and unconstrained cube.
-
Unavailable for bounded cube, LPM2, and SCPS. Their
supported analysis route is
as_svrepdesign(type = "subbootstrap"), a generic PPS bootstrap approximation that does not recreate the constraints or the spatial algorithm.
Exact recursive formulas exist for generalized Brewer (Brewer 2002, ch. 9) but are \(O(N^3)\) and impractical for large frames.
Units with \(\pi_i = 1\) (certainty selections) are handled internally. The joint matrix is computed on the stochastic part and reassembled with \(\pi_{kl} = 1\) for certainty-certainty pairs and \(\pi_{kl} = \pi_l\) for certainty-stochastic pairs.
To use the joint matrix for single-stage variance estimation, pass it
via pps = survey::ppsmat(joint_matrix). See
?joint_expectation for the full method-by-method
breakdown.
Rows and columns follow first appearance in the sample. A WR stage contains one row per distinct selected population unit rather than one row per draw. At a later stage below a WR parent, each parent draw occurrence defines its own conditional block. Cross-occurrence entries are products of marginal conditional chances. Pair the matrix with stage-specific unit identities, not blindly with every descendant sample row.
3. With-replacement row structure
For WR methods (srswr, pps_multinomial),
each draw produces one row with .draw_k indexing the hit.
At export, .draw_k becomes the sampling-unit identifier in
the ids formula, and the FPC is set to \(\infty\). The survey package then uses the
Hansen–Hurwitz variance estimator, which is the standard approach for WR
designs.
4. Chromy / PMR semantics
pps_chromy is treated identically to WR at export (no
FPC, no pps argument). The exact Sen–Yates–Grundy formula
for PMR requires \(E(n_i) E(n_j) - E(n_i
n_j)\) as pairwise weights (Chromy 2009, eq. 5), but
survey::ppsmat() reads \(\pi_i\) from the diagonal, which for PMR
contains \(E(n_i^2)\), not \(E(n_i)\). The Hansen–Hurwitz approximation
is conservative and more stable (Chauvet 2019).
5. Panel assignments
The .panel column is not passed to survey. Panels are
randomized fixed-quota partitions for rotation or workload distribution.
The full-sample weights in .weight are correct for the
combined sample. Panel assignment is not an additional
probability-sampling phase, so a single panel does not automatically
have inclusion probability \(\pi_i /
k\). Conditional on the frozen quotas, its probability is the
block’s realized quota over the block size. Multiplying one panel’s
weights by \(k\) is therefore not
generally valid for population inference.
6. Certainty flag
samplyr stores .certainty_k as a logical
column. At export, certainty units (\(\pi_i =
1\)) are placed in a synthetic stratum
(.cert_stratum = "certainty"), added to the strata formula
alongside user-defined strata. The certainty stratum contributes zero
variance (it is a census) and does not inflate degrees of freedom (Cochran 1977, ch. 11;
Särndal, Swensson, and Wretman 1992,
ch. 3.5).
7. Design metadata
samplyr stores the full sampling_design
object, stage labels, seed, and execution history as attributes of the
tbl_sample. These are accessible via
get_design() and get_stages_executed() but are
not round-tripped through the survey export.
Panel partitioning semantics
There are two distinct operations here and they have different
statistical status. Partitioning
(panels = k with no schedule) splits one execution into
\(k\) groups for rotation or workload.
It does not change any inclusion probability and does not change
.weight. Wave activation (a scheduled
panels plus execute(wave = t)) selects a
subset of those groups, and that is a second sampling phase
with its own inclusion probability and its own weight factor. The rules
below are the partitioning rules, while activation is covered under
“Wave activation” after them.
Neither operation replenishes the sample. A programme that draws
refreshment cohorts against later frame vintages is a registry of
separate executions, not a partition of one. See
rotation_program() and
vignette("rotating-panels").
execute(..., panels = k) partitions the sample into
\(k\) non-overlapping rotation groups
by randomized fixed quota inside frozen ordered blocks:
- Assignment pools: each selection stratum of the assignment stage is a pool, cut into consecutive blocks of \(k \lceil 2 / r_{\min} \rceil\) assignment units, where \(r_{\min}\) is the fewest panels any declared wave activates. Without a schedule, \(r_{\min} = 1\) and the block is \(2k\). A schedule keeping two of four panels live blocks at 4 rather than 8. The take of \(r_{\min}\) panels must leave at least two units per block, which is the smallest take carrying a within-block variance estimate. Every block carries a fixed quota per panel and its labels are permuted within the block, so each unit carries each panel with probability \(1/k\) and panel sizes within a pool differ by at most one.
- Clustered designs: panels are assigned at the assignment stage’s cluster or occurrence level, then propagated to every unit below. By default that stage is the first, so a PSU carries one panel and all its members inherit it. Under a with-replacement stage the assignment unit is the realized draw, so one population cluster drawn twice may carry two panels.
-
Which stage assigns:
panel_stagemoves the assignment to another stage, and every later stage still inherits its ancestor’s panel. Assigning at stage 1 rotates whole primary units. Assigning lower down rotates units inside parents that stay in the survey, which is the address-panel design. -
What moving it costs: a pool becomes the stage’s
own strata inside each realized parent and never crosses a parent, so
pools are smaller and
small_poolbinds more often. Selection certainty then counts only at the assignment stage, in both directions: a certainty primary unit does not make its households permanent, and a certainty selection below the assignment stage does not keep its household in every wave. -
Control sorting: the
controlargument todraw()orders the units before they are cut into blocks, so every panel inherits the same spread over that order. - Certainty units: labelled from their own pools. A certainty unit is in the sample at every occasion and consumes no rotating quota.
The block sizes and the realized block-by-panel quotas are recorded
with the sample and written by write_design(), because
activating a subset of the panels is a simple random subsample without
replacement of the block quota, not of \(1/k\).
Moving the assignment below stage 1 is not a default and should not
be treated as one. Holding the parent fixed while its members rotate can
bias cross-sectional estimates over time, and a unit that cannot move
between parents is balanced for net change but not for gross change.
vignette("rotating-panels") works through the case.
Weights are not adjusted by partitioning. Panels structure the output but do not affect the selection mechanism or inclusion probabilities. Analyze the combined sample with the stored weights.
Wave activation
Activating a subset of the panels is the second operation, and it is
a probability subsample rather than a filter.
execute(master, wave = t) returns a two-phase sample: the
master is phase 1, the activation is phase 2, the blocks are the phase-2
strata and their frozen quotas the phase-2 population counts. Each
weight is multiplied by the inverse of its unit’s activation
probability, which is its block’s frozen quota for the active panels
over the block size. That factor is exact and is generally not
the nominal \(k / r\) implied by
counting groups, because a block that took the remainder carries a
different quota. Certainty units are activated with probability one and
their weights are untouched.
This is why the block sizes and realized block-by-panel quotas are recorded with the sample rather than recomputed: without them the subsample would have to be assumed to be \(1/k\), which it is not.
PRN and sample coordination semantics
samplyr passes permanent random numbers through to
sondage without transformation. PRN-compatible methods
(bernoulli, pps_poisson, pps_sps,
pps_pareto) produce a deterministic selection given the PRN
values. When PRN is supplied, the seed argument to
execute() does not affect selection.
Positive and negative coordination across survey waves is a
user-level workflow, not a package feature. samplyr passes
PRNs to the selection method, the user manages their persistence and
update. See vignette("sampling-coordination"), including
the negative-coordination update \(u_{\text{new}} = (u - \pi) \bmod 1\) (Ohlsson
1995).