Constrained stratified allocation

Distribute a total sample size across strata defined by a single stratification variable, under a fixed total n, target CV, or budget. When the design uses multiple stratification variables (e.g. region and urbanicity), cross them into a single variable beforehand so that each row of frame represents one unique stratum.

Usage

n_alloc(frame, ...)

# Default S3 method
n_alloc(
  frame,
  ...,
  domains = NULL,
  n = NULL,
  cv = NULL,
  budget = NULL,
  measures = NULL,
  targets = NULL,
  psu = NULL,
  certainty_cutoff = NULL,
  n_psu_per_zone = NULL,
  objective = NULL,
  alloc = c("neyman", "optimal", "proportional", "power"),
  unit_cost = NULL,
  alpha = 0.05,
  deff = 1,
  resp_rate = 1,
  df = NULL,
  min_n_stratum = NULL,
  alloc_q = 0.5,
  fpc = c("unit", "stage", "none"),
  plan = NULL
)

# S3 method for class 'svyplan_prec'
n_alloc(frame, ..., n = NULL, cv = NULL, budget = NULL)

Arguments

frame

For the default method: a stratum-level data frame describing the population you want to sample. Each row represents one stratum, a subgroup of the population defined by a stratification variable such as region, age group, or urbanicity. The values in this frame typically come from a census, a population register, or a previous survey. Any stratum table with the columns below works, including a pool summary carried back from an executed sample once the measure columns are added.

When a design stratifies by several variables at once (e.g. region $\times$ urbanicity), cross them into a single variable before calling n_alloc (e.g. with interaction()) so that each row maps to exactly one population cell.

When measures and targets are supplied, frame instead describes atomic allocation strata for joint constrained allocation. It requires unique, non-missing stratum and positive finite N columns. It may contain unit_cost, max_weight, take_all, and complete domain classification columns named by targets$domain. In that mode sd, var, mean and p belong in measures instead. Each row must be homogeneous for every domain classification used by a target. For a fixed-take two-stage design, also supply whole N_psu and n_per_psu, plus positive cost_psu and cost_ssu. A three-stage design additionally requires whole N_ssu and n_per_ssu, and positive cost_tsu. N_psu is the number of PSUs in the stratum. In three-stage mode, N_ssu is the aggregate number of SSUs across all PSUs in the stratum (not the number per PSU). Stage costs replace unit_cost.

Required columns:

N

Number of units (e.g. households, individuals) in each stratum. These are population counts, not sample sizes. Must be positive and finite.

sd or var (legacy mode)

A measure of how spread out the variable of interest is within each stratum. Provide exactly one:

sd: the stratum standard deviation ($\sqrt{\text{variance}}$), or
var: the stratum variance.

Both must be non-negative and finite. When all strata have equal variability (or variability is unknown), a constant column (e.g. sd = 1) yields proportional-to-size allocation.

Optional columns:

stratum: A label identifying each stratum (e.g. "Urban", "Rural"). If omitted, row numbers are used. Must be unique, or unique within each domain when domains is set.
mean or p: The stratum population mean or proportion of the variable of interest. Required when solving for cv or using alloc = "power", because the coefficient of variation is defined relative to the mean. Use mean for continuous variables and p (in $[0, 1]$) for binary (yes/no) variables.
alloc_measure: Positive finite stratum measure used only by Bankier allocation (alloc = "power"). Defaults to N * abs(mean) (or N * p). Population CV is derived from sd / abs(mean) and requires a nonzero mean in every stratum.
unit_cost: Per-unit interviewing cost in each stratum (positive, finite). Set higher values for strata that are more expensive to reach. Defaults to 1 everywhere (equal cost).
max_weight: Maximum allowed sampling weight $N_h / n_h$. Caps how under-represented a stratum can be. Use NA for strata without a cap. It bounds the gross weight, the units drawn, not the response-adjusted $N_h/(n_hr_h)$, which differs once resp_rate varies by stratum.
take_all: Logical (or 0/1). If TRUE, every unit in the stratum is included, a census stratum. Useful for small strata whose total population is tiny enough to enumerate.
deff, resp_rate: Per-stratum design effect and expected response rate, for designs whose fieldwork differs across strata. NA in either column falls back to the scalar argument. The deff and resp_rate arguments override these columns when given as vectors. In joint allocation a non-missing measures row overrides them instead.

For svyplan_prec objects: a precision result from prec_alloc().

...

Additional arguments passed to methods. Unused arguments are rejected.

domains

Character vector of column names in frame to treat as domain identifiers, or NULL (default) for no domains. All names must exist in frame. Domains define sub-populations that each contain one or more strata. When cv is the target, precision is enforced within every domain (see Details).

n

Total sample size. Specify exactly one of n, cv, or budget.

cv

Target coefficient of variation (relative standard error). For example, cv = 0.05 means the standard error of the estimated population mean or total should be at most 5 percent of the estimate. Requires mean or p in frame. When domain columns are present, this target is enforced in each domain. Specify exactly one of n, cv, or budget.

budget

Total field budget. In the classic modes, specify exactly one of n, cv, or budget, and budget caps the affordable allocation. In joint constrained allocation it must be accompanied by objective: the allocation minimizes the priority-weighted relative variance of the objective components subject to the budget, every hard target in targets, and the stratum bounds. Supplying budget without objective in joint mode is an error, because the two readings are incompatible.

measures

Optional long data frame for joint constrained allocation, with one row per required stratum and indicator name. Each row must contain either p in [0, 1], or mean and exactly one of non-negative sd or var. In element-level designs, p implies the finite-population variance N * p * (1 - p) / (N - 1), using its stratum's N > 1, consistently with prec_prop(). Explicit var and sd values are used as supplied. Multistage designs, including PSU-register designs, retain the working proportion variance p * (1 - p). Optional row-specific deff and resp_rate values take precedence, and NA falls back to the stratum's frame column of the same name, then to the scalar argument. Indicator values (p, mean, sd, var) belong here, and frame columns of those names are refused. In fixed-take multistage mode, icc_psu is required and var_ratio_psu defaults to 1. Three-stage mode also requires icc_ssu and derives var_ratio_ssu as var_ratio_psu * (1 - icc_psu) when it is absent, which is the value the variance decomposition implies. Stage parameters may instead be supplied as stratum defaults in frame, in which case non-missing measure-row values take precedence. Only rows selected by a target and its domain are value-validated. Unused rows are ignored after their non-empty, unique keys and frame-stratum membership are checked. Must be supplied with targets.

targets

Optional long data frame with one precision requirement per row. It requires indicator name and exactly one positive cv, moe, or rmoe. Optional domain and level select a frame domain (.overall and level = NA denote the whole population), alpha overrides the scalar value for MOE, and constraint supplies a stable unique identifier. Must be supplied with measures. Optional when objective and budget are given, which requests the best design a budget can buy with no hard precision requirement.

psu

Optional PSU register for certainty-aware allocation, one row per PSU, with stratum and N, the count of ultimate units in that PSU. The PSU sizes must sum to frame$N within each stratum, and every PSU must hold at least its stratum's n_per_psu unless flagged certainty. Optional psu_id is carried through and optional certainty forces a PSU into the certainty part whatever its size. It is an alternative to the frame column N_psu, which says only how many PSUs a stratum has, and requires the frame column n_per_psu and measures$icc_psu. A stratum with a single PSU has no between-PSU variance, so its icc_psu is taken as 0 and may be NA. Joint constrained allocation to targets only, two stages only, and without budget. See the certainty section under Details.

certainty_cutoff

Inclusion probability at which a PSU of the psu register is taken with certainty, in $(0, 1]$. One value applies to every stratum, and one per stratum, as a vector or as a frame column of the same name, relaxes the rule where it is set. The argument overrides the column. NULL (default) uses the column where present and 1 otherwise, which is the rule without a cutoff. Requires psu. See the certainty section under Details.

n_psu_per_zone

NULL (default), 1 or 2, the PSUs drawn from each zone when each stratum's remainder is cut into zones of about equal size. With 2, every stratum that draws from its remainder has an estimable variance. With 1, the variance is estimated by collapsing zones in the groups $psu$.pair fixes before selection. An optional psu$zone_order column sets the order zones are cut in, such as a geographic order, and size, largest first, is the default. Requires psu. See the certainty section under Details.

objective

Optional estimates whose priority-weighted relative variance is minimized among the allocations that meet targets and the budget. Either a character vector of indicator names (overall domain, equal priority) or a data frame with name plus optional domain, level, non-negative priority, and component columns, using the same indicator-domain identifiers as targets. At least one priority must be positive. Priorities need not be normalized. Requires budget and must be supplied with measures.

alloc

Allocation rule: "neyman" (default), "optimal", "proportional", or "power".

unit_cost

Optional scalar or length-nrow(frame) vector of per-stratum unit costs, overriding frame$unit_cost. Not used with fixed-take multistage joint allocation, which requires stage costs.

alpha

Significance level, default 0.05.

deff

Design effect multiplier (> 0). A scalar applies to every stratum. A length-nrow(frame) vector gives one per stratum and overrides a deff column in frame.

resp_rate

Expected response rate, in (0, 1]. Default 1. A scalar applies to every stratum. A length-nrow(frame) vector gives one per stratum and overrides a resp_rate column in frame. See the nonresponse section of svyplan-package for what this adjustment does and does not claim.

df

Degrees of freedom of the variance estimator the allocation will have, typically sampled PSUs minus strata, and available from design_df(). It switches the quantile used to translate a moe target in and to report moe out. A cv target carries no quantile and is unaffected. NULL (default) applies no adjustment.

min_n_stratum

Optional minimum sample size per stratum.

alloc_q

Bankier allocation exponent, used only when alloc = "power". Numeric scalar in $[0, 1]$. With the default alloc_measure, alloc_q = 1 gives Neyman allocation. At alloc_q = 0 weights depend on population CVs and the response/design-effect adjustment, independently of alloc_measure. This does not generally equalize domain CVs. Default 0.5.

fpc

Which finite population correction the variance carries in cluster allocation: "unit" (default) for the ultimate-unit 1 - n / N on the whole inflated variance, "stage" for the stage-by-stage correction, exact under the equal-size, equal-take variance-component model this path assumes, or "none" for n_cluster()'s with-replacement model. "stage" requires N_psu, and N_ssu as well at three stages. Refused outside cluster allocation and with a psu register. See the correction section in Details.

plan

Optional svyplan() object providing design defaults.

Value

A svyplan_n object with type = "alloc" and a stratum-level allocation table in $detail (also available via as.data.frame()), with columns:

stratum, N, sd, unit_cost: Stratum identifiers and inputs carried over from the frame. In cluster mode unit_cost is the derived effective per-element cost cost_psu / n_per_psu + cost_ssu (1 when no stage costs were given), while sd stays the input stratum SD.
n: Allocated sample size (continuous).
n_int: Integer allocation. In n mode an element design keeps the requested total (bounded largest-remainder rounding). A cluster design fields whole PSUs at one whole take per stratum, near its cost-optimal value, so its total can miss n by at most half the largest take. The allocation is re-solved within what each stratum can field at its take, and each rounded PSU count then moves by at most one PSU, jointly across strata, toward n. A take is raised above that value only when the strata at their N_psu limits would otherwise fall more than half a take short. In cv mode each stratum is rounded up so the integer design meets the target. In budget mode units are added by variance reduction per unit cost so the integer design stays within budget. Always inside the integerized bounds (ceiling(.lower), floor(.upper)), except that an N_psu cap bounds the whole PSU count rather than the element count, since the whole take can round up. An error is raised when no integer allocation satisfies them, and when no whole-unit cluster design comes within half a take of n, above or below. That error gives the nearest total the design can field, the smallest when lower bounds and one PSU per stratum force more, the largest when N and N_psu allow less.
weight: Base weight N / n of the issued sample, before nonresponse or calibration adjustment.
n_eff: Effective sample size n * resp_rate / deff.
.lower, .upper: Bounds applied to the stratum (from min_n_stratum, max_weight, take_all, N_psu, or N).
.binding: Whether the allocation sits on one of its bounds.
.bound_source: Which constraint produced the bound the allocation sits on, and NA when it sits on none.
.psu_frac: Cluster mode with N_psu only: the share of the stratum's PSU universe the design selects, n_psu / N_psu.
.se, .moe, .rmoe, .cv: Precision of the stratum's own mean under the allocation, on the same scale as mean (or p). Multiply .se by N for the stratum total. .rmoe is .moe as a fraction of the stratum mean. .cv and .rmoe are NA unless the frame supplies mean or p.
.share: The stratum's share of the design variance of the overall mean. Sums to 1 and identifies where precision is actually being bought. It is NA when the design variance is zero, which a census is, since no stratum then holds a share of nothing.
take_all, mean: Present when take-all strata or stratum means were supplied.
n_per_psu, n_psu: Cluster mode only: the continuous per-stratum take and implied number of PSUs (n / n_per_psu).
n_psu_int, n_per_psu_int: Cluster mode only: the whole-unit field design (whole PSUs and whole takes), chosen so that budget designs stay within budget and cv designs meet the target. Here n_int = n_psu_int * n_per_psu_int.
n_psu, n_psu_int, n_per_psu, n_per_ssu: Joint fixed-take multistage mode: continuous and whole first-stage decisions and the fixed later-stage takes. n_per_ssu is present only for three stages. Public n remains in ultimate-unit units, so n = n_psu * n_per_psu for two stages and n = n_psu * n_per_psu * n_per_ssu for three stages. The same exact identities hold for n_int and n_psu_int.
n_psu_certain, n_psu_rest, threshold: With psu only: the number of certainty PSUs, the number of PSUs left in the remainder, and the certainty threshold certainty_cutoff * n_per_psu / (n / N) at the returned allocation.
n_certain_int, n_psu_draw: With psu only: the units the operational design fields in the certainty PSUs, and the number of remainder PSUs it draws. Here n_int = n_certain_int + n_psu_draw * n_per_psu.
n_zone: With n_psu_per_zone only: the zones the remainder is cut into, so n_psu_draw = n_zone * n_psu_per_zone.

The result also carries an $operational list describing the integer field design. In joint mode it contains total n, variable field cost, the full target table in constraints, all_pass, and the number of integer-repair additions. In legacy modes it also reports the applicable aggregate se, moe, and cv. Top-level quantities describe the continuous design, while as.integer() returns the operational total.

With psu, the result adds $psu, one row per PSU of the register:

psu_id, stratum, N: Carried over from the register. psu_id is present when the register has it.
certainty: Whether the PSU is in the certainty part.
n_take: The take the operational design fields in the PSU. A certainty PSU carries its own whole take at the stratum rate, capped at its size, and a remainder PSU carries n_per_psu.
.certainty_source: Why a certainty PSU is certainty. "threshold" when its size reaches the threshold, "supplied" when psu$certainty flagged it below the threshold, "operational" when it is below the threshold but the whole-PSU remainder draw of $operational would give it an inclusion probability of one, and "orbit" when the resolution of a cycle holds it although no rule requires it. NA for a remainder PSU. summary() counts them.
.threshold, .distance: The stratum's threshold, and the relative distance N / .threshold - 1. A PSU with a small .distance is the one whose classification the allocation can tip.
.zone: With n_psu_per_zone only: the remainder PSU's zone within its stratum, NA for a certainty PSU.
.pair: With n_psu_per_zone = 1 only: the variance group of the PSU's zone, numbered across strata, NA outside any zone.

$optimization$certainty reports how the classification was resolved:

verdict: "converged" when the classification repeated itself, "cycle" when it alternated between classifications, and "limit_reached" when the iteration limit came first.
orbit, iterations: The number of classifications in the repeating orbit (1 when converged), and the passes made before it repeated.
settle_iterations: The solves used to settle the allocation under the held classification.
absorbed: PSUs added to the held classification because the settled allocation, or the repaired field design, made them certainty.
feasible_member: Whether the held member of the orbit met every target before settling. When none did, the most certain member is held and the settle step closes the gap.
fixed_point: Whether the returned classification is exactly the one the returned allocation implies under the threshold and the whole-PSU draw, with any psu$certainty flags. It is FALSE only when PSUs have .certainty_source equal to "orbit", and the printed plan says so.

In fixed-budget objective mode the result adds $objective, one row per component with its priority, relative variance .relvar, equivalent .cv, weighted .contribution, and .share of the total, plus the scalar $objective_value. $operational carries the same two fields for the integer design together with budget_residual. $params$mode is "budget_objective" rather than "targets", and $optimization adds the budget, whether it binds, its residual, the local budget_sensitivity (change in objective value per extra unit of budget), and the epsilon-constraint objective_bound and its multiplier.

Details

Joint constrained allocation

Supplying both measures and targets requests a minimum-cost allocation meeting every indicator-domain precision requirement at once, and adding objective with budget asks instead for the best design that budget can buy. This generalized mode is mutually exclusive with scalar n, cv, budget, domains and an explicit alloc rule, because the target table defines the constraints and the allocation is solved jointly. Its input contract, its multistage variance model, its infeasibility messages and the predict() frontier are on n_alloc-generalized.

For fieldwork, use $detail$n_int, not the analytical continuous n. Fixed-take count columns are numeric vectors whose values are exact whole numbers. The _int suffix denotes operational integer semantics rather than R's 32-bit integer storage type. Result lists contain diagnostics and are not promised as a long-term serialization format.

Building the frame

The frame is a data frame where each row is one stratum of your target population. It summarizes what you know about each subgroup before sampling. A typical workflow:

Identify strata from a census or register (e.g. provinces, urban/rural areas, age groups).
Look up N: the population count per stratum.
Estimate sd: the standard deviation of your key variable within each stratum (from a pilot survey, a previous census, or expert judgement). If unknown, set sd = 1 everywhere for proportional allocation.
Add mean or p if you want to solve for a target CV.

A minimal frame:

frame <- data.frame(
  stratum = c("Urban", "Rural"),
  N       = c(50000, 120000),
  sd      = c(12, 20)
)

When a design stratifies by several variables (e.g. region $\times$ urbanicity), cross them into one variable first:

frame$stratum <- interaction(frame$region, frame$urban, drop = TRUE)

This ensures that each row maps to exactly one population cell and that the allocation formulas apply to the correct per-stratum N and sd pairs.

Cluster designs within strata

Adding a icc_psu column to the frame turns the allocation into a stratified two-stage design (e.g. enumeration areas then households within each stratum). Under the cluster variance model the problem reduces to the element allocation above with the stratum SD inflated to sd * sqrt(var_ratio_psu * (1 + icc_psu * (n_per_psu * resp_rate - 1))) and, when stage costs are given, a per-element cost of cost_psu / n_per_psu + cost_ssu. The clustering penalty is paid on the take that responds, n_per_psu * resp_rate, since a unit that does not respond contributes no within-cluster observation. The cost is paid on the gross take, since it is issued either way. The two coincide at resp_rate = 1. The whole-cluster operational search reads the same responding take, so a design it accepts is one the continuous reduction also accepts. All solve modes, allocation methods, and constraints work unchanged. The n, cv, and budget modes keep their meanings.

Cluster-mode columns:

icc_psu (required): within-PSU homogeneity per stratum, e.g. from varcomp() with strata. When PSUs are selected with probability proportional to size, pass prob to varcomp() as well.
var_ratio_psu (optional, default 1): variance ratio per stratum.
n_per_psu (optional): fixes the per-stratum take, in the operational whole-unit design as much as in the continuous one. Any NA entries are replaced by the cost-optimal take sqrt(cost_psu / cost_ssu * (1 - icc_psu) / icc_psu). Holding the take is what n mode gives up its exact total for. The stratum is fielded as a whole number of clusters of that size, so its element count lands on the nearest multiple.
cost_psu, cost_ssu (together): per-PSU and per-element costs. These are required for budget mode or when n_per_psu is not fixed. They replace unit_cost, which is not allowed in this mode.
N_psu (optional): the number of PSUs available in the stratum.

What `N_psu` does, and what it does not

N_psu is a feasibility constraint and nothing more. It caps the allocation at N_psu * n_per_psu ultimate units and caps the whole-unit design at N_psu clusters, matching the bound the fixed-take path already applies. Leaving it out preserves the unbounded behavior, in which the allocation may ask for more PSUs than a stratum contains.

On its own it does not activate a first-stage finite population correction. Under the default fpc = "unit" no factor keyed to the PSU sampling fraction enters the variance, however large a share of the PSU universe the design takes, and print() says so once that share is appreciable. Supplying N_psu bounds the design. Asking for the correction as well is what fpc = "stage" does.

So at n_psu == N_psu the default still carries between-PSU variance even though every PSU has been selected. That is deliberately conservative and is no longer a literal variance representation. For the same reason, a target this path reports as unreachable at the PSU bound may be reachable under a finite-population first stage, and the error says so rather than claiming the precision is impossible. Under fpc = "stage" the between-PSU component does reach zero there, and the whole variance reaches zero only when every stage is enumerated.

Because taking every PSU leaves the within-PSU take in force, it does not enumerate a stratum, and take_all is refused in cluster mode for the same reason the fixed-take path refuses it.

The whole-unit design keeps its PSU count within N_psu as well. When the strata at their N_psu limits cannot hold n at their cost-optimal whole takes, a take is raised above its continuous value. When no whole take brings the design within half a take of n, n_alloc() stops with an error giving the largest total it can field.

Which correction the variance carries

fpc chooses, and the three choices are nested rather than arbitrary. Write $f_1$ for the PSU sampling fraction and $f_2$ for the within-PSU one, both measured on what responds. Under the equal-size, equal-take variance-component model this path assumes, the two-stage without-replacement variance corrects each component by its own stage:

$$V = (1-f_1)\frac{S_1^2}{a} + (1-f_2)\frac{S_2^2}{am}.$$

"unit" (default) applies the ultimate-unit 1 - n / N to the whole clustering-inflated variance, so the between-PSU component is attenuated by it too. Since $n/N = f_1f_2$, that factor is at least as large as either exact one, so this overstates both components rather than either being understated. It is conservative against the exact without-replacement variance at every sampling fraction.
"stage" applies the expression above, and its three-stage analogue. It needs N_psu, and N_ssu as well at three stages. It is exact for that model and not beyond it: unequal PSU sizes, a size-proportional first stage, and stochastic rather than expected nonresponse each need a more general variance expression, so it is refused with a psu register.
"none" applies no correction at any stage, which is n_cluster()'s model exactly. Use it to compare the two interfaces at a sampling fraction where the correction would otherwise separate them.

So none $\ge$ unit $\ge$ stage at every sampling fraction, and the three coincide as the fractions vanish.

The default stays "unit" deliberately. "stage" is exact under that model, but it buys the exactness by depending on the mean PSU size $M = N/N_{psu}$, which the with-replacement form does not use at all. Under unequal PSU sizes that mean is a stand-in, and planners usually know it less well than they know icc_psu and the take. Choose "stage" when the PSU sizes really are close to equal and the fractions are large enough to matter.

Note that the two fractions are controlled by different quantities. $f_1$ falls as the design takes a smaller share of the PSU universe, but $f_2$ is the realized take over $M$ and does not move with N at all. A large population reached through a few large PSUs leaves $f_2$ wherever it was, so "stage" is not in general a small correction to "none".

fpc is refused outside cluster allocation, where there are no stages to choose between and the correction is always 1 - n / N.

Because icc_psu already accounts for the clustering, leave deff at 1 unless it captures a different source of design effect (e.g. weighting loss). A clustering deff on top of icc_psu would double-count. The constraints min_n_stratum and max_weight stay in element units. For fielding, use the whole-unit design in n_psu_int and n_per_psu_int (n_int = n_psu_int * n_per_psu_int). Its actual field cost and precision are reported in $operational and, in budget mode, never exceed the budget.

Certainty PSUs from a register

N_psu gives only the number of PSUs in a stratum. Supply psu to give their sizes instead. It has one row per PSU with stratum and N, and the sizes must sum to frame$N in every stratum. This two-stage joint-allocation mode requires n_per_psu in frame and icc_psu in measures.

A remainder PSU is drawn for the take, so a PSU holding fewer units than n_per_psu is refused, named by its psu_id or its row. Merge it with a neighbouring PSU before planning, as is usual when PSUs are formed (Valliant, Dever and Kreuter 2018, sec. 10.3) and as merge_psus() does, or flag it in certainty, where it is subsampled at the stratum's sampling fraction. A psu_id must name each PSU once.

icc_psu describes the remainder, whose PSUs are selected with probability proportional to size. It is the size-weighted share of the stratum variance that lies between PSUs, which varcomp() returns when prob is each PSU's share of its stratum's units, as in the last example. A stratum with a single PSU has no variance between PSUs, so its icc_psu is taken as 0 whatever is supplied and may be NA. Every other stratum needs a value in $[0, 1]$.

With take $b_h$, sampling fraction $f_h = n_h / N_h$ and certainty_cutoff $c_h$, a PSU is certain when its inclusion probability reaches $c_h$, that is when its size reaches

$$N_{hi} \ge c_h b_h / f_h.$$

The remainder is fielded as a whole number of PSUs drawn with probability proportional to size, and rounding that number up raises every remainder PSU's inclusion probability. A PSU whose probability reaches $c_h$ under that draw is certainty too, and the draw is recomputed without it until no PSU reaches $c_h$, the iterative rule for selection proportional to size. A PSU below the threshold can be certainty this way, and so can a stratum's only PSU. A remainder always draws at least one PSU, even when the rounded certainty takes use up the stratum's allocation, so no part of a stratum is left out of the sample. The operational design is therefore one any fixed-size proportional-to-size method can draw as it stands.

The default $c_h = 1$ takes with certainty only the PSUs the design cannot sample less often. A cutoff below one, such as the 0.80 Valliant, Dever and Kreuter (2018, sec. 9.6) mention, takes more PSUs with certainty and removes their clustering from the design effect, so the sample usually falls. The classification is discrete, so where the loop cycles the plan moves in steps as the cutoff changes and need not fall with it. Price the choice with predict.svyplan over certainty_cutoff.

n_psu_per_zone = 2 cuts each stratum's remainder, sorted by psu$zone_order or by size, into zones of about equal total size and draws two PSUs from each, the deep stratification of area samples (Valliant, Dever and Kreuter 2018, sec. 9.5 and Example 3.13). The remainder draw rounds up to a whole number of zones, which can add one PSU in a stratum. The rule above then applies inside each zone. A PSU is certain when its probability within its zone reaches $c_h$, and so is every PSU of a zone, or of a remainder, that holds no more PSUs than it draws. A stratum that draws from its remainder then draws at least two PSUs, so its variance can be estimated. The design effect does not credit the zones, so the planned precision stays conservative.

n_psu_per_zone = 1 draws one PSU from each zone, the finest stratification of the remainder. A zone holding a single PSU is a census, so more PSUs are certainty and the sample tends to be the smallest of the three settings. No zone yields two PSUs, so the variance needs zones collapsed, and $psu$.pair fixes the groups before selection. Zones are paired in zone order within a stratum, with three together where the count is odd. Strata whose remainder is a single zone are grouped with each other in frame row order, so order the rows to put similar strata next to each other. One draw per zone has the with-replacement variance, so the groups are analysed without a finite population correction. In that form the collapsed estimator overestimates the variance, by more as the zone totals in a group differ (Valliant, Dever and Kreuter 2018, sec. 15.5.3). A without-replacement correction such as Brewer's can turn it into an underestimate when the zones in a group are alike.

The certainty part has no first-stage sampling variance. The remainder has the usual clustering component, and the two are combined into the anticipated design effect. psu$certainty can add a PSU to the certainty part. It cannot remove a PSU either rule makes certainty.

The continuous allocation uses an aggregate planning approximation. The noncertainty part uses a with-replacement clustering approximation, a standard approach in PPS sample-size planning. Each part receives its proportional share of the stratum allocation, so with $N_h^C$ units in certainty PSUs, $N_h^R$ in the remainder, $\delta_h$ for icc_psu and response rate $r_h$, the design effect is

$$d_h = \frac{N_h^C + N_h^R \{1 + \delta_h (b_h r_h - 1)\}}{N_h}.$$

It depends on which PSUs are certainty, not on the allocation. A common stratum-level finite population correction is then applied. The clustering reads the responding take $b_h r_h$, as the N_psu form does, while the threshold and the draw read the issued take $b_h$. A PSU-level response rate and a var_ratio_psu other than 1 have no place in this design effect and are refused.

The threshold and allocation determine each other, so the solver iterates. $optimization$certainty records whether it converged, cycled, or reached its iteration limit. $detail gives the split and threshold by stratum, while $psu gives each PSU's classification, its source, and n_take, the take the operational design fields in it. A certainty PSU carries its own whole take at the stratum rate, capped at its size, and a remainder PSU carries n_per_psu. The takes are the operational design's own numbers, so sum(n_take[certainty]) = n_certain_int in every stratum.

A fit solved with a register keeps it in $params$psu. Its presence is the stable test for a certainty-aware fit, and it's possible to key on it rather than on $detail column names.

The operational design is fieldable: certainty PSUs use their whole take at the stratum rate and the remainder uses whole PSUs at n_per_psu. Thus n_int = n_certain_int + n_psu_draw * n_per_psu in every stratum. $operational reports its count, cost and precision. The rounded certainty takes $q_{hi}$ run above the stratum rate and the remainder gets what is left, so the field design is assessed on its own takes and draw $a_h$, each part a stratum of its own:

$$V_h = S_h^2 \Big[\sum_{i \in C} N_{hi}^2 \Big(\frac{1}{r_h q_{hi}} - \frac{1}{N_{hi}}\Big) + (N_h^R)^2 \{1 + \delta_h (b_h r_h - 1)\} \Big(\frac{1}{r_h a_h b_h} - \frac{1}{N_h^R}\Big)\Big].$$

At the stratum rate the expression equals the aggregate one above. Where the field design misses a target, whole remainder PSUs (whole zones when zoned) are added one at a time, each where it most reduces the failing targets' excess variance per unit of cost, as the solver's own integer repair does. Strata outside the failing targets' domains are left alone. If the remainder feeding a failing target is used up, or no addition reduces it, n_alloc() stops with an error naming each failing target with its achieved and required values. That error concerns this repair under the held classification and certainty takes, not every possible field design.

Use cost_psu and cost_ssu together to price PSU visits and interviews. A certainty PSU's visit is paid whatever its take, so the solver prices an interview at cost_ssu + (N_h^R / N_h) * cost_psu / n_per_psu, and $params$achieved$cost and $optimization$cost add the certainty visits as a fixed cost, as do the stratum costs summary() reports. The optimality certificate, its gap and its multipliers refer to the marginal-cost problem without that constant, which has the same optimum, and hold for the classification the plan settled on, not across classifications. A budget is refused with a register, because the field design rounds its draws up and does not hold certainty visits within a budget. predict.svyplan can compare fixed values of n_per_psu.

Domains vs. strata

Domains are specified via the domains parameter. Domain columns partition strata into sub-populations. Each domain groups one or more strata. When cv is specified, the algorithm finds the minimum total $n$ such that the worst-case domain CV meets the target, i.e. every domain achieves the required precision.

In n or budget mode, domains affect reporting only: per-domain precision metrics appear in $domains but the allocation itself treats all strata globally.

Allocation methods

Allocation is controlled by the alloc parameter (same methods as strata_bound()):

proportional: $n_h \propto N_h / r_h$
neyman: $n_h \propto N_h S_h \sqrt{d_h / r_h}$
optimal: $n_h \propto N_h S_h \sqrt{d_h / r_h} / \sqrt{c_h}$
power: Bankier (1988) allocation, $n_h \propto (S_h/|\mu_h|) X_h^q\sqrt{d_h/r_h}$, where $X_h$ is frame$alloc_measure, defaulting to $N_h|\mu_h|$. Every stratum needs a finite nonzero mean or p. The population CV $S_h/|\mu_h|$ is distinct from the target cv, which describes the relative standard error of an estimate. Zero-variance strata have zero allocation weight and remain subject to the allocation bounds. If all variances are zero, the existing fallback allocates in proportion to population size, with a warning. Use explicit domain targets when comparable domain precision is required.

The design effect $d_h$ and response rate $r_h$ enter only when they vary by stratum. A value shared by every stratum is a constant factor and cancels out of a proportional weighting, leaving the classical rules above it. Neyman and optimal allocation carry $\sqrt{d_h/r_h}$ from minimizing $\sum W_h^2S_h^2d_h/(r_hn_h)$ under a constraint on issued count or cost, respectively. The power rule retains that adjustment while weighting population CVs by alloc_measure^alloc_q. Proportional carries $1/r_h$ instead, and no $d_h$ at all, because it is a count rule, not a variance optimum, and its purpose is a self-weighting sample, so it is the responding sample it holds proportional to $N_h$.

Stratum allocations are rounded to integers using the ORIC method (Cont and Heidari, 2015). Constraints (min_n_stratum, max_weight, take_all) are enforced via recursive Neyman allocation (RNA, Wesolowski et al., 2021).

When budget is specified, the algorithm finds the maximum affordable allocation under unit costs.

References

Valliant, R., Dever, J. A., & Kreuter, F. (2018). Practical Tools for Designing and Weighting Survey Samples (2nd ed.). Springer. Chapters 5 and 9.

Bankier, M. D. (1988). Power allocations: determining sample sizes for subnational areas. The American Statistician, 42(3), 174–177.

Bethel, J. (1989). Sample allocation in multivariate surveys. Survey Methodology, 15(1), 47–57.

Cont, R. and Heidari, M. (2015). Optimal rounding under integer constraints. arXiv preprint arXiv:1501.00014.

Wesolowski, J., Wieczorkowski, R. and Wojciak, W. (2021). Optimality of the recursive Neyman allocation. Journal of Survey Statistics and Methodology, 10(5), 1263–1275.

Examples

frame <- data.frame(
  stratum = c("A", "B", "C"),
  N    = c(4000, 3000, 3000),
  sd   = c(10, 15, 8),
  mean = c(50, 60, 55),
  unit_cost = c(1, 1.5, 1)
)

n_alloc(frame, n = 600)
#> Stratum allocation (neyman, 3 strata)
#> field design: n = 600, cv = 0.0079, cost = 724
#> continuous optimum: n = 600, cv = 0.0079, se = 0.4305
#> (deff = 1)
#> design df = 597
n_alloc(frame, cv = 0.03)
#> Stratum allocation (neyman, 3 strata)
#> field design: n = 46, cv = 0.0294, cost = 56
#> continuous optimum: n = 44.23479, cv = 0.0300, se = 1.6350
#> (deff = 1)
#> design df = 43

frame_constraints <- transform(
  frame,
  max_weight = c(25, 20, NA),
  take_all = c(FALSE, FALSE, TRUE)
)

n_alloc(frame_constraints, budget = 3500, alloc = "optimal", min_n_stratum = 40)
#> Stratum allocation (optimal, 3 strata)
#> field design: n = 3403, cv = 0.0076, cost = 3500
#> continuous optimum: n = 3403.425, cv = 0.0076, se = 0.4125
#> (min_n_stratum = 40, deff = 1)
#> design df = 401

frame_domains <- data.frame(
  province = c("North", "North", "South", "South"),
  stratum = c("Urban", "Rural", "Urban", "Rural"),
  N    = c(2000, 3000, 1800, 3200),
  sd   = c(12, 18, 10, 16),
  mean = c(55, 48, 58, 50)
)

n_alloc(frame_domains, domains = "province",
       cv = 0.04, alloc = "power", alloc_q = 0.3)
#> Stratum allocation (power, 4 strata)
#> field design: n = 110, cv = 0.0271, cost = 110
#> continuous optimum: n = 108.2371, cv = 0.0273, se = 1.4160
#> (deff = 1)
#> design df = 106
#> Domains: 2
#> 
#>  province .domain .n       .se     .moe     .rmoe      .cv    .cost
#>  North    5_North 58.67031 2.03200 3.982647 0.07839856 0.0400 59   
#>  South    5_South 49.56679 1.97265 3.866324 0.07311505 0.0373 50   

# Joint targets for two indicators and an overlapping domain
joint_frame <- data.frame(
  stratum = c("North urban", "North rural", "South urban"),
  region = c("North", "North", "South"),
  N = c(1000, 1800, 1200)
)
joint_measures <- data.frame(
  stratum = rep(joint_frame$stratum, 2),
  name = rep(c("coverage", "income"), each = 3),
  p = c(0.5, 0.4, 0.6, NA, NA, NA),
  mean = c(NA, NA, NA, 50, 55, 60),
  sd = c(NA, NA, NA, 10, 12, 15)
)
joint_targets <- data.frame(
  name = c("coverage", "income"),
  domain = c(".overall", "region"),
  level = c(NA, "North"),
  cv = c(0.06, NA),
  moe = c(NA, 3)
)
joint <- n_alloc(
  joint_frame, measures = joint_measures, targets = joint_targets
)
prec_alloc(joint, n = joint$detail$n_int)
#> Joint allocation precision (2 constraints)
#> targets: all pass

# The best design a fixed budget can buy, subject to the same targets
bought <- n_alloc(
  joint_frame, measures = joint_measures, targets = joint_targets,
  objective = "income", budget = 4000
)
bought$objective
#>         component   name   domain level priority .relvar .cv .contribution
#> 1 income@.overall income .overall  <NA>        1       0   0             0
#>   .share
#> 1     NA

# Unequal priorities, and an objective on a specific domain
priorities <- data.frame(
  name     = c("income", "coverage"),
  domain   = c(".overall", "region"),
  level    = c(NA, "North"),
  priority = c(2, 1)
)
n_alloc(
  joint_frame, measures = joint_measures, targets = joint_targets,
  objective = priorities, budget = 4000
)
#> Joint constrained allocation (Bethel)
#> question: best design affordable within a budget of 4000
#> field design: n = 4000, cost = 4000 (2 targets, all pass)
#> continuous optimum: n = 4000, cost = 4000
#> objective: weighted relative variance 7.63819e-20 continuous, 7.63819e-20 operational
#> the budget is not binding: the allocation sits at its upper bounds
#> 
#>  component             priority .cv          .share
#>  income@.overall       2        0.000000e+00 0     
#>  coverage@region=North 1        2.763728e-10 1     
#> 
#> active bounds: 0 lower, 3 upper

# What would other budgets buy? (cost-versus-objective frontier)
predict(bought, data.frame(budget = c(3000, 4000, 5000)))
#>   budget    n cost objective_value          cv n_int cost_int .binding
#> 1   3000 3000 3000    3.907646e-06 0.001976777  3000     3000     TRUE
#> 2   4000 4000 4000    0.000000e+00 0.000000000  4000     4000    FALSE
#> 3   5000 4000 4000    0.000000e+00 0.000000000  4000     4000    FALSE
#>   .feasible
#> 1      TRUE
#> 2      TRUE
#> 3      TRUE

# The same target tables with fixed two-stage takes
joint_frame$N_psu <- c(100, 150, 100)
joint_frame$n_per_psu <- c(8, 10, 12)
joint_frame$cost_psu <- c(300, 400, 450)
joint_frame$cost_ssu <- c(25, 30, 35)
joint_measures$icc_psu <- rep(c(0.03, 0.05, 0.08), 2)
joint_measures$var_ratio_psu <- 1
joint_cluster <- n_alloc(
  joint_frame, measures = joint_measures, targets = joint_targets
)
joint_cluster$detail[, c("stratum", "n_psu_int", "n_per_psu", "n_int")]
#>       stratum n_psu_int n_per_psu n_int
#> 1 North urban        11         8    88
#> 2 North rural        17        10   170
#> 3 South urban        11        12   132

# Fixed three-stage takes: PSUs, SSUs within PSU, ultimate units within SSU
joint_frame$N_ssu <- c(600, 1000, 750) # aggregate SSUs by stratum
joint_frame$n_per_ssu <- c(3, 4, 3)
joint_frame$cost_tsu <- c(5, 6, 7)
joint_measures$icc_ssu <- rep(c(0.10, 0.08, 0.12), 2)
joint_measures$var_ratio_ssu <- 1
joint_three <- n_alloc(
  joint_frame, measures = joint_measures, targets = joint_targets,
  min_n_stratum = 20
)
joint_three$detail[, c(
  "stratum", "n_psu_int", "n_per_psu", "n_per_ssu", "n_int"
)]
#>       stratum n_psu_int n_per_psu n_per_ssu n_int
#> 1 North urban         7         8         3   168
#> 2 North rural         9        10         4   360
#> 3 South urban         6        12         3   216
prec_alloc(joint_three, n = joint_three$detail$n_int)
#> Joint allocation precision (2 constraints)
#> targets: all pass

# Stratified two-stage design (EAs then households per stratum)
frame_cluster <- data.frame(
  stratum   = c("Urban", "Rural"),
  N         = c(50000, 150000),
  sd        = c(0.45, 0.48),
  mean      = c(0.35, 0.25),
  icc_psu = c(0.03, 0.08),
  cost_psu  = c(300, 600),
  cost_ssu  = c(40, 60)
)

n_alloc(frame_cluster, cv = 0.05)
#> Stratum allocation (neyman, two-stage, 2 strata)
#> field design: n = 2021, n_psu = 171, cv = 0.0498, cost = 206500
#> continuous optimum: n = 1979.882, cv = 0.0500, se = 0.0138
#> (deff = 1)
#> design df = 169

# Certainty PSUs from a register that lists every PSU with its size.
# The "City" stratum is a single PSU.
set.seed(1)
size <- c(round(rlnorm(60, log(150), 0.8)), 4000)
register <- data.frame(
  psu_id = seq_along(size),
  stratum = rep(c("North", "South", "City"), c(30, 30, 1)),
  N = size
)
pop <- data.frame(
  stratum = rep(register$stratum, size),
  psu_id = rep(register$psu_id, size),
  y = rbinom(sum(size), 1, rep(plogis(rnorm(61, -0.5, 0.4)), size))
)
# icc_psu under size-proportional selection: each PSU's share of its stratum
pop$share <- ave(pop$y, pop$psu_id, FUN = length) /
  ave(pop$y, pop$stratum, FUN = length)
vc <- varcomp(y ~ psu_id, strata = ~stratum, prob = ~share,
              data = subset(pop, stratum != "City"))

reg_frame <- aggregate(N ~ stratum, register, sum)
reg_frame$n_per_psu <- 10
reg_measures <- aggregate(cbind(p = y) ~ stratum, pop, mean)
reg_measures$name <- "y"
# NA for "City": a single PSU has no between-PSU variance
reg_measures$icc_psu <- vc$strata$icc_psu[
  match(reg_measures$stratum, vc$strata$stratum)
]
reg_measures
#>   stratum         p name    icc_psu
#> 1    City 0.3330000    y         NA
#> 2   North 0.3906302    y 0.03284821
#> 3   South 0.3829961    y 0.03343549

reg_fit <- n_alloc(
  reg_frame, measures = reg_measures,
  targets = data.frame(name = "y", cv = 0.05), psu = register
)
reg_fit$detail[, c("stratum", "n_psu_certain", "n_psu_draw", "n_int")]
#>   stratum n_psu_certain n_psu_draw n_int
#> 1    City             1          0   162
#> 2   North            12          9   264
#> 3   South            10         11   274

# Relaxing the certainty rule to an inclusion probability below one
predict(reg_fit, data.frame(certainty_cutoff = c(1, 0.9, 0.8)))
#>   certainty_cutoff n_psu_certain n_psu_draw        n     cost n_int cost_int
#> 1              1.0            23         20 696.2411 696.2411   700      700
#> 2              0.9            27         16 686.9169 686.9169   692      692
#> 3              0.8            30         13 680.1932 680.1932   686      686
#>   .feasible
#> 1      TRUE
#> 2      TRUE
#> 3      TRUE

# No zones, one PSU per zone, and two PSUs per zone
predict(reg_fit, data.frame(n_psu_per_zone = c(NA, 1, 2)))
#>   n_psu_per_zone n_psu_certain n_psu_draw n_zone        n     cost n_int
#> 1             NA            23         20      0 696.2411 696.2411   700
#> 2              1            39          7      7 665.1471 665.1471   674
#> 3              2            30         14      7 680.1932 680.1932   696
#>   cost_int .feasible
#> 1      700      TRUE
#> 2      674      TRUE
#> 3      696      TRUE