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,
  seed = 1L
)

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.

The frame is checked against the digest recorded at execution. Other columns, another row order or other column types are accepted when every pool keeps the size and selection chances the sample was drawn with. Rows added or removed, a changed design value or a missing design column are refused with samplyr_error_joint_frame_mismatch, and so is any difference at all for a stage whose method depends on row order (systematic, pps_systematic, pps_chromy, registered methods), whose matrix is that of the order. Under frame_digest = "none" there is nothing to compare against, and a warning says so.

...

These dots are for future extensions and must be empty. stages and the arguments after it 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, nsim and seed, 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. Analytic methods ignore it. Raising it narrows the simulation error, at a cost linear in nsim.

seed

Single integer seeding the simulated methods (default 1), which analytic methods ignore. The session's random stream is restored afterwards, so the result is an exact function of the sample, nsim and seed, and calling this never moves a later execute(). At the default nsim a chromy matrix carries simulation error of a few percent, so two seeds give slightly different answers, neither more correct than the other.

Value

With waves, a tibble with one row per block of the frozen assignment, described in the "Activation mode" section. 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 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 are the stage's sampled units in the order they first appear in the sample. For a sample with one row per stage-1 unit, the stage-1 matrix therefore follows the sample rows, which is how survey::ppsmat() reads it. 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.

Each matrix is dense and covers every sampled unit of the stage, across pools as well as within them, so m sampled units take $8 m^2$ bytes: about 80 MB at 3,200 units, 0.8 GB at 10,000 and 12.8 GB at 40,000. For a large stage, request only the stages needed with stages.

Details

Without frame, the computation uses the frame digest recorded at execution, so a sample that traveled without its (possibly confidential) frame still yields joint expectations at the method's stated quality. This needs each pool's exact chances, which the digest always holds for cluster stages and constant-chance element stages, but for element stages with varying chances only under execute(frame_digest = "full"). A summarized digest refuses rather than approximates, so pass the frame instead. An exact record of approximate targets does not make their probabilities exact. A supplied frame must keep the pools the sample was drawn from, as frame describes, and its units must be uniquely identifiable within each stratum/cluster group by their column values.

A sample drawn from separately supplied stage registers needs either its intact digest or all of those registers, in the order they were supplied. One frame is refused with samplyr_error_frame_count, because a lower-stage register holds no rows for the population an upper stage selected from.

Each PPS stage's full-population first-order quantities are rebuilt from its method and measure of size, with per-stratum targets (n_h) replayed from the allocation execute() used (proportional, Neyman, optimal, etc.), so they match those at sampling time. The matching sondage joint function then gives the submatrix of sampled units.

For stratified or conditional (within-cluster) stages, joint quantities are computed independently within each group, and cross-group entries are products of the marginal chances. The rows are then put in the order the units first appear in the sample, so the units of one group need not be contiguous. Below a WR parent, pair the matrix with stage-specific identities in this order, not blindly with every sample row when descendants duplicate a selected unit.

With certainty selections ($\pi_i = 1$) in a WOR design, the joint probabilities of the 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.

Exact vs. approximate computation

The accuracy of the returned matrix depends on the sampling method, as the two tables below show. variance_estimation covers how a matrix enters a variance estimate.

WOR methods ($\pi_{kl}$)

Every matrix here comes from sondage::joint_inclusion_prob().

Method	Quality
`pps_cps`	Exact, Aires' formula
`pps_sampford`	Exact, Sampford design
`pps_systematic`	Exact, circular-interval overlap
`pps_poisson`	Exact, independent draws
`pps_brewer`	Approximate$^*$, high-entropy
`pps_sps`	Approximate, high-entropy
`pps_pareto`	Approximate, high-entropy
`cube`	Approximate when unconstrained
`lpm2`	Unavailable, spatial spreading not represented
`scps`	Unavailable, spatial spreading not represented

Poisson draws are independent, so $\pi_{kl} = \pi_k \pi_l$. Unconstrained cube uses the high-entropy approximation too.

Systematic PPS commonly has zero pair probabilities. Its exact matrix describes the design, but a zero pair probability rules out a design-unbiased variance estimator, and near-zero pairs make it unstable.

$^*$ Exact recursive formulas for Brewer's joint inclusion probabilities (Brewer 2002, ch. 9) are $O(N^3)$, impractical beyond a few hundred units, whereas the high-entropy approximation is $O(N^2)$. It assumes the design is close to the maximum-entropy design with the same marginal $\pi_i$ (Hajek 1964; Brewer and Donadio 2003), which strong ordering, balancing or extreme probabilities can break, so validate variance and interval coverage for the intended design and population. SPS and Pareto also enter it with approximate first-order targets.

Bounded cube, LPM2, and SCPS designs are refused, because count constraints and spatial spreading alter pairwise selection beyond the approximation. variance_estimation gives their replicate variance.

A stage that draws a certainty plan's zones treats each zone as a draw of its own, so PSUs of different zones have the product of their chances. It needs frame, since the digest does not record the zone mates of a selected PSU. A plan drawing one PSU per zone is refused (samplyr_error_joint_method_unsupported). Two PSUs of one zone are never drawn together, so no unbiased variance estimator can use its matrix.

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

Both matrices come from sondage::joint_expected_hits().

Method	Quality
`pps_multinomial`	Exact, analytic
`pps_chromy`	Approximate, Monte Carlo, 10 000 draws

For pps_multinomial, $E(n_k \cdot n_l) = n(n-1) p_k p_l + n p_k \mathbf{1}_{k=l}$. The sequential dependence of pps_chromy admits no closed form for $E(n_k \cdot n_l)$, so sondage estimates it by Monte Carlo simulation under nsim and seed.

Activation mode

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. joint_expectation(master, waves = c(t, s)) reads it from the record, with neither the frame nor a simulation.

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 has one row per block and stays small at any sample size, yet every entry of the full matrix is recoverable from it, since within a block the joints take only two values. The conditional covariance kernel is block-diagonal, not the joint-probability matrix: across blocks the joint is $p_i p_j$, generally not zero, while the covariance is zero.

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

These are joint expectations of the activation indicators, conditional on the phase-1 units and the frozen quotas, not the unconditional joint inclusion probabilities of the 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 neither of the two waves needs to have been materialized.

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 conditional Poisson (CPS) joint probabilities: Aires, N. (1999). Algorithms to find exact inclusion probabilities for conditional Poisson sampling and Pareto $\pi$ps sampling designs. Methodology and Computing in Applied Probability, 1(4), 457-469. doi:10.1023/A:1010091628740

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

# A single-stage stratified Sampford sample, whose joint matrix is exact
sample <- sampling_design() |>
  stratify_by(region) |>
  draw(n = 5, method = "pps_sampford", mos = households) |>
  execute(bfa_eas, seed = 2025)

jip <- joint_expectation(sample, bfa_eas)

# survey reads the matrix by row, one row per sampled unit
svy <- as_svydesign(sample, pps = survey::ppsmat(jip[[1]]))
survey::svytotal(~population, svy)
#>               total     SE
#> population 20854071 297547

# 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>