Compute pairwise joint expectations from a sample

Reconstructs the second-order design quantities for PPS stages. For without-replacement (WOR) stages, this produces the joint inclusion probabilities $\pi_{kl}$. For with-replacement (WR) and PMR stages, this produces the joint expected hits $E(n_k \cdot n_l)$.

Usage

joint_expectation(
  x,
  frame = NULL,
  ...,
  stages = NULL,
  waves = NULL,
  nsim = 10000L
)

Arguments

x: A tbl_sample object produced by execute().
frame: The frame originally passed to execute(). A data frame is the one shared frame. An ordered list of data frames is the stage registers, one per executed stage, in the order they were supplied to execute(). Either way each frame must contain the columns its stage sampled on (strata variables, cluster variables, measure of size). When NULL (the default), the computation uses the frame digest recorded on the sample instead.
...: These dots are for future extensions and must be empty. stages and the arguments after it follow ..., so each must be named exactly: the singular stage is reported rather than prefix-matched.
stages: An integer vector of stage numbers to compute, or NULL (default) to compute all PPS stages. Non-PPS stages produce NULL entries in the returned list.
waves: Two wave numbers of a scheduled master, which selects activation mode: the result describes how those two occasions of the rotation overlap rather than how the master was selected. Give the same wave twice for the joint expectation within one wave. NULL (default) keeps the stage behavior. Mutually exclusive with stages, frame and nsim, none of which activation mode uses.
nsim: Positive integer number of simulations used for Chromy's pairwise expected hits (default 10000). Also forwarded to registered WR joint_fns that explicitly declare an nsim formal. Ignored by analytic methods.

Value

With waves, a tibble with one row per block of the frozen assignment. See "Activation mode" below. Otherwise a named list of length equal to the number of executed stages. Each element is either:

For PPS WOR stages: a square matrix of joint inclusion probabilities $\pi_{kl}$, usable with survey::ppsmat() for exact variance estimation.
For PPS WR/PMR stages (pps_multinomial, pps_chromy): a square matrix of joint expected hits $E(n_k \cdot n_l)$.
NULL for non-PPS stages (SRS, systematic) or stages not requested via the stages argument.

Rows and columns represent stage-specific sampled units in first appearance order. At a WR stage, repeated hits of the same population unit appear once, so dimensions match the number of distinct sampled units (or clusters). At a later stage below a WR parent, each parent draw occurrence defines a separate conditional block, so the same child population identity can appear in more than one block.

Details

Without frame, the computation runs off the frame digest recorded at execution: the digest holds each pool's exact resolved chance vector, which is all the joint computation needs, so a sample that traveled without its (possibly confidential) frame still yields exact joint expectations. This requires an exact chance representation: cluster stages always have one, element stages with constant chances have one, and element stages with varying chances keep one only under execute(frame_digest = "full"). A summarized representation refuses rather than approximates. With frame, the quantities are replayed against it as before. The frame must be unchanged since execution (validate_frame() reports drift).

A sample drawn from separately supplied stage registers needs either its intact digest or the ordered list of those registers. One lower-stage register cannot reconstruct upper-stage joint quantities: it holds no rows for the population the upper stage selected from.

For each PPS stage, the function:

Reconstructs the full-population first-order quantities from the frame using the stage's method and measure of size
Dispatches to the appropriate sondage joint probability or joint expected hits function
Extracts the submatrix corresponding to sampled units

For stratified stages, the target sample size per stratum (n_h) is reconstructed by replaying the same allocation logic used during execute() (proportional, Neyman, optimal, etc.) against the frame. This ensures first-order quantities match what was computed at sampling time, regardless of allocation method.

For stratified or conditional (within-cluster) stages, joint quantities are computed independently within each group. Blocks follow their first appearance in the sample, as do units within a block. Cross-block entries are products of the corresponding marginal chances.

A stage below a WR parent is conditional on each parent draw occurrence, not only on the parent's population identity. Repeated hits of one parent therefore produce separate independent child blocks. Pair the returned matrix with stage-specific identities in this order. Do not pair it blindly with every sample row when descendants duplicate a selected unit.

Exact vs. approximate computation

The accuracy of the returned matrix depends on the sampling method. Some algorithms yield closed-form joint probabilities. Others require approximation or simulation.

WOR methods ($\pi_{kl}$)

samplyr method	sondage function	Quality
`pps_cps`	`joint_inclusion_prob()`	Exact (Aires' formula via C)
`pps_sampford`	`joint_inclusion_prob()`	Exact (Sampford design)
`pps_systematic`	`joint_inclusion_prob()`	Exact (circular-interval overlap)
`pps_poisson`	`joint_inclusion_prob()`	Exact ($\pi_{kl} = \pi_k \pi_l$, independent draws)
`pps_brewer`	`joint_inclusion_prob()`	Approximate$^*$ (high-entropy / Hajek-Brewer-Donadio)
`pps_sps`	`joint_inclusion_prob()`	Approximate (high-entropy / Hajek-Brewer-Donadio)
`pps_pareto`	`joint_inclusion_prob()`	Approximate (high-entropy / Hajek-Brewer-Donadio)
`cube`	`joint_inclusion_prob()`	Approximate when unconstrained (high-entropy / Hajek-Brewer-Donadio)
`lpm2`	unavailable	Spatial spreading is not represented
`scps`	unavailable	Spatial spreading is not represented

$^*$ Exact recursive formulas for Brewer's joint inclusion probabilities exist (Brewer 2002, ch. 9) but are $O(N^3)$, making them impractical for frames of more than a few hundred units. The high-entropy approximation is $O(N^2)$ and sufficiently accurate for variance estimation in practice. The same trade-off applies to SPS and Pareto, whose exact joint probabilities would require enumerating the combinatorial sample space.

The high-entropy approximation assumes the design is close to the maximum-entropy design with the same marginal $\pi_i$ (Hajek 1964; Brewer and Donadio 2003). This is a good approximation for most PPS designs and is the same quantity that underlies the Berger (2004) variance estimator used by survey::svydesign(pps = "brewer"). For CPS (conditional Poisson / maximum entropy), the joint probabilities are exact by definition.

Bounded cube, LPM2, and SCPS designs are rejected because count constraints and spatial spreading alter pairwise selection behavior beyond the available approximation. Use as_svrepdesign() with type = "subbootstrap" for a generic bootstrap approximation instead.

WR/PMR methods ($E(n_k \cdot n_l)$)

samplyr method	sondage function	Quality
`pps_multinomial`	`joint_expected_hits()`	Exact (analytic: $n(n-1) p_k p_l + n p_k \mathbf{1}_{k=l}$)
`pps_chromy`	`joint_expected_hits()`	Approximate (Monte Carlo simulation, 10 000 replicates)

For pps_chromy, the sequential dependence structure does not admit a closed-form expression for $E(n_k \cdot n_l)$. sondage uses Monte Carlo simulation (default 10 000 replicates) to estimate the pairwise expectations. Increasing nsim reduces Monte Carlo error at the cost of computation time.

Activation mode

joint_expectation(master, waves = c(t, s)) answers a different question from the stage modes above. A scheduled master assigned every unit to a panel at its draw and froze the block-by-panel quotas, so how two occasions of the rotation overlap is already determined and needs neither the frame nor a simulation. It is read from the record.

Conditional on the frozen quotas, the panels inside a block of m assignment units are an arrangement of that block's labels. Writing $a_t$ and $a_s$ for the units the panels active at each wave take from the block, and $a_\cap$ for those active at both:

$$P(i \in W_t) = a_t / m$$ $$P(i \in W_t, i \in W_s) = a_\cap / m$$ $$P(i \in W_t, j \in W_s) = \frac{a_t a_s - a_\cap}{m (m - 1)}$$

for units $i \neq j$ of one block, and the product of the marginals for units of different blocks, whose arrangements are drawn independently. The within-wave case is the third expression at t == s.

The result is one row per block, and stays small at any sample size. Within a block the joints take only two values, and across blocks they factorize as products of the marginals, so every entry of the full matrix is recoverable from this table. It is the conditional covariance kernel, not the joint-probability matrix, that is block-diagonal: across blocks the joint is $p_i p_j$, which is generally not zero, while the covariance is.

pool, stratum, class and block identify the block. units is $m$. take_1, take_2 and take_both are the three takes. prob_1 and prob_2 are the marginals. joint_same and joint_distinct are the two joint expectations. has_pair is FALSE for a block of one unit, where no distinct pair exists and joint_distinct is NA rather than zero.

Quotas are frequently unequal, because a pool that is not a multiple of the block size gives one block an extra unit, so these are read from the record rather than derived from the panel count. A certainty block is permanent and takes every unit at every wave, which makes joint_distinct exactly one.

What this states, and what it does not. These are joint expectations of the activation indicators, conditional on the phase-1 units and on the frozen quotas. They are not the unconditional joint inclusion probabilities of the complete two-phase design, which also carry the master's own pairwise term. How the two combine for a variance of change is not settled here. as_svydesign() carries the activation as a second phase for ordinary totals.

Activation mode takes a master, not a materialized wave, so a pair of waves neither of which has been materialized can be asked for.

Limitations

The frame-free path requires a digest with exact chances: the default summary digest suffices for cluster stages and constant-chance element stages. Element stages with varying chances need execute(frame_digest = "full"). Otherwise pass the frame.
When frame is supplied it must be unchanged from what was passed to execute(), and units in it must be uniquely identifiable within each stratum/cluster group by their column values.
A sample drawn from one register per stage needs all of them, in the order they were supplied. One frame is refused with samplyr_error_frame_count rather than computed from the wrong population.
For WOR designs with certainty selections ($\pi_i = 1$), the joint matrix is decomposed: certainty units are separated from the stochastic part, the joint probabilities for non-certainty units are computed from the reduced $\pi$ vector, and the full matrix is reassembled with $\pi_{ij} = 1$ for certainty pairs and $\pi_{ij} = \pi_j$ for certainty x non-certainty pairs.

References

High-entropy approximation: Hájek, J. (1964). Asymptotic theory of rejective sampling with varying probabilities from a finite population. Annals of Mathematical Statistics, 35(4), 1491-1523.

Brewer, K.R.W. and Donadio, M.E. (2003). The high entropy variance of the Horvitz-Thompson estimator. Survey Methodology, 29(2), 189-196.

Exact Brewer joint probabilities: Brewer, K.R.W. (2002). Combined Survey Sampling Inference: Weighing Basu's Elephants. Arnold, ch. 9.

The variance estimator these underlie: Berger, Y.G. (2004). A simple variance estimator for unequal probability sampling without replacement. Journal of Applied Statistics, 31(3), 305-315.

Examples

sample <- sampling_design() |>
  add_stage() |>
    stratify_by(region) |>
    cluster_by(ea_id) |>
    draw(n = 5, method = "pps_brewer", mos = households) |>
  add_stage() |>
    draw(n = 12) |>
  execute(bfa_eas, seed = 2025)
#> Warning: Stage 2: selected every unit available in the pools it executed.
#> ✖ Requested 780 units, selected all 65 available.
#> ℹ Exhausted pools: "33131", "33250", "10550", "36826", "44761", and 60 more.
#> ℹ This stage contributes no sampling variance. Earlier stages are unaffected:
#>   the design as a whole is a census only if every stage is.

# Compute joint probabilities for stage 1
jip <- joint_expectation(sample, bfa_eas, stages = 1)

# Use with survey package for exact variance (WOR stages)
svy <- as_svydesign(sample, pps = survey::ppsmat(jip[[1]]))
#> Warning: Exact PPS variance (`pps`) is single-stage in survey.
#> ℹ Exporting the stage-1 design only; later-stage sampling variance is not
#>   represented.
#> ℹ Omit `pps` for multi-stage linearization with Brewer's approximation at the
#>   PPS stage.

# Compute all PPS stages at once
jip_all <- joint_expectation(sample, bfa_eas)

# A sample drawn from one register per stage passes them as a list
regions <- dplyr::distinct(bfa_eas, region, .keep_all = TRUE)
registers <- sampling_design() |>
  add_stage() |>
    cluster_by(region) |>
    draw(n = 3, method = "pps_brewer", mos = households) |>
  add_stage() |>
    draw(n = 12) |>
  execute(regions, bfa_eas, seed = 2025)

jip_registers <- joint_expectation(registers, list(regions, bfa_eas))
# Activation mode: how far two occasions of a rotation overlap. Four
# panels rotate two at a time, so each is live for two consecutive waves.
rotation <- data.frame(
  panel = rep(1:4, times = 4),
  wave = rep(1:4, each = 4),
  active = c(
    TRUE, TRUE, FALSE, FALSE,
    FALSE, TRUE, TRUE, FALSE,
    FALSE, FALSE, TRUE, TRUE,
    TRUE, FALSE, FALSE, TRUE
  )
)

master <- sampling_design() |>
  draw(n = 40) |>
  execute(bfa_eas, seed = 2025, panels = rotation)

# Waves 1 and 2 share one panel of the two each activates.
joint_expectation(master, waves = c(1, 2))
#> # A tibble: 10 × 16
#>     pool stratum class    activation block wave_1 wave_2 units take_1 take_2
#>    <int> <chr>   <chr>    <chr>      <int>  <int>  <int> <int>  <int>  <int>
#>  1     1 NA      rotating rotating       1      1      2     4      2      2
#>  2     1 NA      rotating rotating       2      1      2     4      2      2
#>  3     1 NA      rotating rotating       3      1      2     4      2      2
#>  4     1 NA      rotating rotating       4      1      2     4      2      2
#>  5     1 NA      rotating rotating       5      1      2     4      2      2
#>  6     1 NA      rotating rotating       6      1      2     4      2      2
#>  7     1 NA      rotating rotating       7      1      2     4      2      2
#>  8     1 NA      rotating rotating       8      1      2     4      2      2
#>  9     1 NA      rotating rotating       9      1      2     4      2      2
#> 10     1 NA      rotating rotating      10      1      2     4      2      2
#> # ℹ 6 more variables: take_both <int>, prob_1 <dbl>, prob_2 <dbl>,
#> #   joint_same <dbl>, joint_distinct <dbl>, has_pair <lgl>

# The same wave twice gives the joint expectation within one wave.
joint_expectation(master, waves = c(2, 2))
#> # A tibble: 10 × 16
#>     pool stratum class    activation block wave_1 wave_2 units take_1 take_2
#>    <int> <chr>   <chr>    <chr>      <int>  <int>  <int> <int>  <int>  <int>
#>  1     1 NA      rotating rotating       1      2      2     4      2      2
#>  2     1 NA      rotating rotating       2      2      2     4      2      2
#>  3     1 NA      rotating rotating       3      2      2     4      2      2
#>  4     1 NA      rotating rotating       4      2      2     4      2      2
#>  5     1 NA      rotating rotating       5      2      2     4      2      2
#>  6     1 NA      rotating rotating       6      2      2     4      2      2
#>  7     1 NA      rotating rotating       7      2      2     4      2      2
#>  8     1 NA      rotating rotating       8      2      2     4      2      2
#>  9     1 NA      rotating rotating       9      2      2     4      2      2
#> 10     1 NA      rotating rotating      10      2      2     4      2      2
#> # ℹ 6 more variables: take_both <int>, prob_1 <dbl>, prob_2 <dbl>,
#> #   joint_same <dbl>, joint_distinct <dbl>, has_pair <lgl>