Skip to contents

The package supports two- and three-stage cluster designs. A primary sampling unit (PSU) is selected at the first stage, a secondary sampling unit (SSU) at the second, and an ultimate unit is the person, household, or other unit on which the outcome is measured. See Getting started for the core sizing and precision workflow and sample allocation for stratified designs.

Multistage cluster designs

For two-stage or three-stage designs, the question is not just “how many?” but “how many clusters and how many units per cluster?” The answer depends on the cost structure and the degree of clustering.

Variance components from frame data

If frame data with cluster identifiers are available, varcomp() estimates between-cluster and within-cluster variance components:

set.seed(1)
n_clust <- 80
n_hh <- 15
cluster_id <- rep(seq_len(n_clust), each = n_hh)
cluster_mean <- rnorm(n_clust, mean = 50, sd = 10)
expenditure <- cluster_mean[cluster_id] + rnorm(n_clust * n_hh, sd = 20)
frame <- data.frame(cluster = cluster_id, expenditure = expenditure)

vc <- varcomp(expenditure ~ cluster, data = frame)
vc
#> Variance components (2-stage)
#> varb = 0.0398, varw = 0.1700
#> icc = 0.1896
#> var_ratio = 1.0589
#> Unit relvariance = 0.1981
as.data.frame(vc)
#>   stages       varb      varw       icc var_ratio unit_relvar
#> 1      2 0.03976226 0.1699763 0.1895801  1.058885   0.1980749

For this two-stage design, icc is the measure of homogeneity \delta = V_b/(V_b + V_w) in the sampling variance decomposition described by Valliant, Dever, and Kreuter (2018). It feeds directly into n_cluster() and the cluster-planning mode of design_effect(). The associated design-effect approximation is k[1 + (b - 1)\delta], where var_ratio supplies k and b is the number of units sampled per PSU.

This measure is related to conventional intraclass correlations and Kish’s rate of homogeneity (roh), but the values are not generally identical. When supplying an ICC from another source, check its definition and the associated design-effect formula. The component ratio above lies between 0 and 1 for nonnegative components with a positive sum. Values near either boundary are special cases for the closed-form optimizer, so n_cluster() rejects values numerically too close to either boundary.

Optimal two-stage allocation

Given per-stage costs and the variance structure, n_cluster() finds the allocation that either minimizes cost for a target CV, or minimizes CV for a given budget:

# Minimize cost to achieve CV = 0.05
n_cluster(stage_cost = c(500, 50), icc = vc, cv = 0.05)
#> Optimal 2-stage allocation
#> field design: n_psu = 26 | n_per_psu = 7 -> total n = 182
#> cv = 0.0496, cost = 22100
#> continuous optimum: n_psu = 26.30386 | n_per_psu = 6.538208 (cv = 0.0500, cost = 21751)
#> design df = 25

# Maximize precision within a budget
n_cluster(stage_cost = c(500, 50), icc = vc, budget = 100000)
#> Optimal 2-stage allocation
#> field design: n_psu = 125 | n_per_psu = 6 -> total n = 750
#> cv = 0.0233, cost = 100000
#> continuous optimum: n_psu = 120.9321 | n_per_psu = 6.538208 (cv = 0.0233, cost = 100000)
#> design df = 124

Passing icc = vc (the varcomp object) automatically extracts icc, unit_relvar, and var_ratio.

Three-stage designs

Add a third cost element and a second icc for three-stage designs:

n_cluster(
  stage_cost  = c(1000, 200, 20),
  icc = c(0.01, 0.05),
  cv    = 0.05
)
#> Optimal 3-stage allocation
#> field design: n_psu = 13 | n_per_psu = 5 | n_per_ssu = 15 -> total n = 975
#> cv = 0.0500, cost = 45500
#> continuous optimum: n_psu = 13.46593 | n_per_psu = 4.974937 | n_per_ssu = 13.78405 (cv = 0.0500, cost = 45333)
#> design df = 12

Response rate in cluster designs

A cluster design can lose units at more than one stage, and the two losses are not interchangeable, so they are named for the stage they act on. resp_rate_psu is the share of selected clusters that can be worked at all, as when whole PSUs are lost to insecurity, flooding, or an inaccessible road. Losing a PSU removes its entire contribution, so it deflates every variance term alike and the first-stage sample is inflated by 1 / resp_rate_psu to cover it.

In budget mode, reducing the PSU response rate alone leaves the allocation unchanged and worsens the achieved CV. Later-stage response rates can change the optimal takes even with a fixed budget. In CV mode, compensate for lost PSUs by increasing the first-stage issue:

n_cluster(stage_cost = c(500, 50), icc = vc, cv = 0.05, resp_rate_psu = 0.90)
#> Optimal 2-stage allocation
#> field design: n_psu = 29 | n_per_psu = 7 -> total n = 203 (net: 183)
#> (resp_rate_psu = 0.90)
#> cv = 0.0495, cost = 24650
#> continuous optimum: n_psu = 29.22651 | n_per_psu = 6.538208 (cv = 0.0500, cost = 24168)
#> design df = 28

Specify response separately for each sampled stage:

Argument Affected unit Interpretation
resp_rate_psu PSU, in either design Probability that a selected cluster can be worked.
resp_rate_ssu SSU, in a three-stage design Response at the second stage, conditional on the PSU being available.
resp_rate Ultimate unit Completion at the last stage, conditional on the earlier stages being available.

For a two-stage design, the SSU is the ultimate unit, so its response belongs in resp_rate. Supplying resp_rate_ssu for that design returns an error.

The three are not interchangeable. Losing whole PSUs deflates every variance term alike, whereas losing ultimate units shrinks the realized cluster and so changes the clustering penalty too. That is why the later-stage rates move the cost-optimal take while the PSU rate does not:

# a DHS-shaped design: EAs, households within EAs, one person per household
n_cluster(stage_cost = c(500, 100, 50), icc = c(0.02, 0.05), cv = 0.05,
          resp_rate_psu = 0.95, resp_rate_ssu = 0.90, resp_rate = 0.85)
#> Optimal 3-stage allocation
#> field design: n_psu = 36 | n_per_psu = 4 | n_per_ssu = 6 -> total n = 864 (net: 628)
#> (resp_rate_psu = 0.95, resp_rate_ssu = 0.90, resp_rate = 0.85)
#> cv = 0.0497, cost = 75600
#> continuous optimum: n_psu = 35.40747 | n_per_psu = 3.689324 | n_per_ssu = 6.686246 (cv = 0.0500, cost = 74438)
#> design df = 35

Sizes and costs stay gross. $n counts the units to issue and $cost pays for them, while the variance reads what they realize. Lumping the later rates into resp_rate_psu overstates the variance, so the plan meets the target at a higher cost than needed and with the wrong take, and the excess cost grows with the between-PSU homogeneity.

Fixed overhead costs

The standard cost model is purely linear: C = c1*n_psu + c2*n_psu*n_per_psu. Real surveys also have a fixed overhead for training, infrastructure, or setup. The fixed_cost parameter adds this term. In budget mode, only budget - fixed_cost is available for variable costs. This reduces the number of PSUs while leaving the continuous cost-optimal take per PSU unchanged. In CV mode, the required allocation stays the same and fixed_cost is added to total cost:

# Budget mode: fixed overhead reduces the allocatable budget
n_cluster(stage_cost = c(500, 50), icc = vc, budget = 100000, fixed_cost = 5000)
#> Optimal 2-stage allocation
#> field design: n_psu = 111 | n_per_psu = 7 -> total n = 777
#> cv = 0.0240, cost = 99350 (fixed: 5000)
#> continuous optimum: n_psu = 114.8855 | n_per_psu = 6.538208 (cv = 0.0239, cost = 100000)
#> design df = 110

# CV mode: same allocation, cost increases by fixed_cost
n_cluster(stage_cost = c(500, 50), icc = vc, cv = 0.05, fixed_cost = 5000)
#> Optimal 2-stage allocation
#> field design: n_psu = 26 | n_per_psu = 7 -> total n = 182
#> cv = 0.0496, cost = 27100 (fixed: 5000)
#> continuous optimum: n_psu = 26.30386 | n_per_psu = 6.538208 (cv = 0.0500, cost = 26751)
#> design df = 25

The same parameter is available when n_cluster() sizes a table of indicators.

Evaluating cluster allocations

Compare a candidate fieldwork allocation with the continuous optimum, keeping the same variance components and costs. Suppose operations can field 100 PSUs with 10 ultimate units each:

alloc <- n_cluster(stage_cost = c(500, 50), icc = vc, budget = 100000)
candidate <- c(n_psu = 100, n_per_psu = 10)
candidate_cost <- candidate[1] * (500 + 50 * candidate[2])

data.frame(
  design = c("Continuous optimum", "100 PSUs with 10 units each"),
  n_psu = c(alloc$n[1], candidate[1]),
  n_per_psu = c(alloc$n[2], candidate[2]),
  cost = c(alloc$cost, candidate_cost),
  cv = c(prec_cluster(alloc)$cv, prec_cluster(n = candidate, icc = vc)$cv)
)
#>                            design    n_psu n_per_psu  cost         cv
#>                Continuous optimum 120.9321  6.538208 1e+05 0.02331895
#> n_psu 100 PSUs with 10 units each 100.0000 10.000000 1e+05 0.02382433

The whole-unit candidate spends the same 100,000 budget. Its CV shows the precision cost of that operational choice relative to the continuous optimum. This comparison assumes that at least 100 PSUs and 10 units per selected PSU are available. Frame limits must also be checked before selecting records.

Vary an uncertain cost before committing to that plan:

predict(alloc, data.frame(cost_psu = c(400, 500, 600)))
#>   cost_psu    n_psu n_per_psu  total_n         cv  cost
#> 1      400 144.4257  5.847951 844.5944 0.02183037 1e+05
#> 2      500 120.9321  6.538208 790.6791 0.02331895 1e+05
#> 3      600 104.3719  7.162248 747.5373 0.02466473 1e+05

The budget stays fixed while the preferred cluster counts and takes change. A higher cost of opening a PSU favors taking more units in each visited PSU.

Vary the budget to see what each increment buys:

predict(alloc, data.frame(budget = c(50000, 100000, 150000, 200000)))
#>   budget     n_psu n_per_psu   total_n         cv   cost
#> 1  50000  60.46604  6.538208  395.3396 0.03297797  50000
#> 2 100000 120.93209  6.538208  790.6791 0.02331895 100000
#> 3 150000 181.39813  6.538208 1186.0187 0.01903984 150000
#> 4 200000 241.86417  6.538208 1581.3583 0.01648899 200000

The cost-optimal take does not depend on the budget, so every increment goes to more PSUs. Without a finite population correction, the CV falls with the square root of the budget, and doubling the budget cuts it by a factor of about 1.4.

Stratified two-stage designs

When a survey stratifies first and clusters within each stratum, use n_alloc() to choose stratum-specific cluster counts and takes. Estimate the variance components within strata from a pilot or previous survey:

set.seed(3)
listing <- data.frame(
  region = rep(c("North", "South"), each = 600),
  ea = rep(1:60, each = 20),
  income = rnorm(1200, rep(c(50, 70), each = 600), 15) +
    rep(rnorm(60, 0, 6), each = 20)
)
vc_region <- varcomp(income ~ ea, data = listing, strata = ~region)
vc_region
#> Variance components (2-stage, 2 strata)
#> 
#>  stratum sd      mean    icc_psu var_ratio_psu varb   varw   unit_relvar
#>  North   15.6736 51.1316 0.1056  1.0489        0.0104 0.0882 0.0940     
#>  South   16.0116 69.3691 0.1595  1.0479        0.0089 0.0469 0.0533

Here ea identifies enumeration areas, the PSUs. Join the components to a planning frame with one row per stratum. N counts ultimate population units, while cost_psu is paid per selected PSU and cost_ssu per ultimate unit:

frame_2stage <- merge(
  data.frame(stratum = c("North", "South"), N = c(40000, 60000)),
  as.data.frame(vc_region)[, c(
    "stratum", "sd", "mean", "icc_psu", "var_ratio_psu"
  )],
  by = "stratum"
)
frame_2stage$cost_psu <- c(400, 550)
frame_2stage$cost_ssu <- c(45, 60)

stratified_plan <- n_alloc(frame_2stage, cv = 0.02)
stratified_plan$detail[, c(
  "stratum", "n_per_psu", "n_int", "n_psu_int", "n_per_psu_int"
)]
#>   stratum n_per_psu n_int n_psu_int n_per_psu_int
#> 1   North  8.675795   135        15             9
#> 2   South  6.951041   203        29             7
stratified_plan$operational[c("n", "cost", "cv")]
#> $n
#> [1] 338
#> 
#> $cost
#> [1] 40205
#> 
#> $cv
#> [1] 0.01966919

n_per_psu is the continuous cost-optimal take. The field design selects n_psu_int whole PSUs at the whole take n_per_psu_int, which gives the n_int ultimate units, and $operational reports that design’s cost and precision. The function chooses cost-optimal takes unless they are fixed by a frame column. take_all is unsupported in this mode, because visiting all PSUs while subsampling within them does not enumerate the stratum.

The take is part of the design, so evaluating the field design, or a variation of it, needs both stage sizes. With them, prec_alloc() reproduces the operational precision:

field <- stratified_plan$detail
prec_alloc(
  stratified_plan,
  n = field$n_int,
  n_per_psu = field$n_per_psu_int
)$cv
#> [1] 0.01966919

Bound the number of available PSUs

Suppose the frame contains only 40 PSUs in the North and 900 in the South:

frame_bounded <- frame_2stage
frame_bounded$N_psu <- c(40, 900)
bounded_plan <- n_alloc(frame_bounded, n = 3000)
bounded_plan$detail[, c(
  "stratum", "n", "n_int", "n_psu_int", "n_per_psu_int", "N_psu",
  ".bound_source"
)]
#>   stratum         n n_int n_psu_int n_per_psu_int N_psu .bound_source
#> 1   North  347.0318   342        38             9    40         N_psu
#> 2   South 2652.9682  2660       380             7   900          <NA>
bounded_plan$operational[c("n", "cost", "cv")]
#> $n
#> [1] 3002
#> 
#> $cost
#> [1] 399190
#> 
#> $cv
#> [1] 0.008653086

The North allocation reaches its PSU ceiling, and .bound_source names the constraint. The remaining sample comes from the South. In whole units the North design stays within its 40 PSUs, and the field total stays within one take of the requested 3,000.

On its own, N_psu is an availability bound. Under the default fpc = "unit" it adds no first-stage finite population correction, so precision keeps the with-replacement first-stage approximation even when a stratum selects most of its PSUs. fpc = "stage" applies the stage-by-stage correction, which needs N_psu, and N_ssu as well at three stages:

c(
  unit = bounded_plan$cv,
  stage = n_alloc(frame_bounded, n = 3000, fpc = "stage")$cv
)
#>        unit       stage 
#> 0.008542285 0.006364018

An unattainable precision target returns an error rather than requesting unavailable PSUs.

Assess the design

Design effects

design_effect() builds the design effect you expect a plan to produce, from the design features you are choosing. The clustering component uses the same variance model as n_cluster() and prec_cluster(), so the two always agree:

design_effect(icc = 0.05, n_per_psu = 20)
#> Planning design effect: 1.9500

Components are selected by the arguments you supply and multiply together. Unequal weighting comes from planned weights or from a planned allocation’s stratum sizes and takes, and the stratification gain (the only component that can fall below 1) from the stratum standard deviations and means:

frame <- data.frame(
  N = c(50000, 120000), n = c(600, 400), sd = c(12, 20), mean = c(55, 48)
)
deff <- design_effect(icc = 0.05, n_per_psu = 20, strata = frame)
deff
#> Planning design effect: 2.6279

The result is a plain number wherever one is expected, and effective_n() is its mirror:

n_prop(p = 0.3, moe = 0.05, deff = deff)
#> Sample size for proportion (wald)
#> n = 848 (p = 0.30, moe = 0.050, deff = 2.63)
#> expected cases = 254.4
effective_n(deff, n = 1000)
#> [1] 380.5354

You can also read the features off a plan you already built, including an n_alloc() allocation:

design_effect(n_cluster(stage_cost = c(500, 50), icc = 0.05, cv = 0.05))
#> Planning design effect: 1.6392

These are planning quantities. After data collection, compute the design effect the sample actually achieved with survey::svymean(..., deff = TRUE), which uses the realized weights, strata, and clusters.

Planned degrees of freedom

For a design-based variance estimator, degrees of freedom depend on the number of independently sampled PSUs and the number of strata. In an unstratified cluster design with m sampled PSUs, the planned count is m - 1:

cluster_plan <- n_cluster(
  stage_cost = c(500, 50),
  icc = 0.05,
  budget = 100000
)
design_df(cluster_plan)
#> Design degrees of freedom (planning)
#> 
#>   df = 79   (80 PSUs - 1 stratum)

design_df() uses the whole-unit PSU count that will be fielded. Pass its result to interval calculations that accept df, such as the Korn–Graubard proportion method. In a stratified design, the count is sampled PSUs minus strata.

This is a planned count. Certainty PSUs do not contribute between-PSU variance, and a domain that cuts across strata needs the realized number of PSUs containing domain members. A plan without that frame information cannot derive those adjustments exactly. See ?design_df for limits by design class.

Multi-indicator multistage

For surveys with several indicators in a cluster design, hand n_cluster() a table of indicators carrying icc_psu columns, along with stage costs:

targets_ms <- data.frame(
  name   = c("stunting", "vaccination"),
  p      = c(0.25, 0.70),
  cv     = c(0.10, 0.08),
  icc_psu = c(0.05, 0.02)
)

n_cluster(indicators = targets_ms, stage_cost = c(500, 50), budget = 80000)
#> Multi-indicator optimal allocation (2-stage)
#> field design: n_psu = 64 | n_per_psu = 15 -> total n = 960
#> worst cv = 0.0729, cost = 80000 (binding: stunting)
#> continuous optimum: n_psu = 67.27204 | n_per_psu = 13.78403 (cv = 0.0728, cost = 80000)
#> 
#>  name        .n        .cv_target .cv_achieved .binding
#>  stunting    491.76043 0.10       0.0728       *       
#>  vaccination  84.08575 0.08       0.0241

For several indicators and domains allocated across explicit strata, continue with Joint allocation. The variance and cost models used here follow Valliant, Dever, and Kreuter (2018).

References

Valliant, Richard, Jill A. Dever, and Frauke Kreuter. 2018. Practical Tools for Designing and Weighting Survey Samples. 2nd ed. Cham: Springer. https://doi.org/10.1007/978-3-319-93632-1.