Strata boundaries for survey design

Determine where to cut a continuous stratification variable to form useful strata. Supports four methods: cumulative root frequency (Dalenius-Hodges), geometric progression, coordinate optimization inspired by Lavallée-Hidiroglou, and Kozak's random search.

Usage

strata_bound(
  x,
  n_strata,
  ...,
  n = NULL,
  cv = NULL,
  method = c("lh", "cumrootf", "geo", "kozak"),
  alloc = c("neyman", "optimal", "proportional", "power"),
  alloc_q = 0.5,
  unit_cost = NULL,
  take_all_above = NULL,
  deff = 1,
  resp_rate = 1,
  n_class = NULL,
  max_iter = NULL,
  n_restart = NULL,
  plan = NULL
)

Arguments

x: Numeric vector: finite stratification variable values. Must not contain missing values.
n_strata: Required integer: number of strata (including take-all if take_all_above is specified). Must be >= 2. Every sampled stratum must hold at least two population units, so a request that only sparse or tied data could satisfy is refused rather than met with a stratum that carries no within-stratum variance. A take-all stratum is enumerated instead of sampled, so one unit is enough there.
...: Unused. Present so that every optional argument must be named. Unused arguments are rejected.
n: Target total sample size, supplied as a whole number. Specify at most one of n or cv. Required for methods "lh" and "kozak".
cv: Target coefficient of variation (relative standard error). For example, cv = 0.05 means the standard error of the estimated population total or mean should be at most 5 percent of the estimate. Specify at most one of n or cv. Required for methods "lh" and "kozak".
method: Stratification method: "cumrootf" (Dalenius-Hodges), "geo" (geometric), "lh" (LH-inspired coordinate optimization), or "kozak" (Kozak's random search). "lh" is a coordinate search in the spirit of the published method, not its recurrence. Default "lh".
alloc: Allocation rule: "proportional", "neyman", "optimal", or "power" (Bankier allocation). Default "neyman". See Details.
alloc_q: Bankier allocation exponent, used only when alloc = "power". Numeric scalar in $[0, 1]$. At alloc_q = 1 the allocation equals Neyman. At alloc_q = 0 weights are proportional to population CVs. This does not generally equalize stratum CVs. Default 0.5.
unit_cost: Per-stratum unit costs, ordered from lowest to highest stratum. Scalar (equal costs) or vector of length n_strata. Default NULL (equal unit costs, $c_h = 1$ for all strata), in which case "optimal" and "neyman" coincide.
take_all_above: Take-all threshold. Units with x >= take_all_above form a census stratum.
deff: Design effect multiplier (> 0), default 1. Inflates the variance of the stratified mean exactly as it does in n_alloc(), so that cv and n mean the same thing in both. A scalar value leaves the boundaries unchanged. See Details.
resp_rate: Expected response rate, in (0, 1], default 1. With a response rate below 1, n is the sample fielded and n * resp_rate the number expected to respond, matching n_alloc().
n_class: Positive whole number of histogram bins. Applies to method = "cumrootf" only. Supplying it with another method is an error rather than a silent no-op. Default NULL (Freedman-Diaconis rule).
max_iter: Positive whole number, the iteration cap. For method = "lh" it counts sweeps over the boundaries (default 200). For "kozak" it caps the moves in each run (default 10,000), which normally stops earlier on its own convergence rule. Supplying it with another method is an error rather than a silent no-op.
n_restart: Positive whole number of Kozak runs from each starting boundary set. Applies to method = "kozak" only. Supplying it with another method is an error rather than a silent no-op. Default NULL (= 5).
plan: Optional svyplan() object providing design defaults. It can supply alloc, alloc_q, deff, resp_rate, and a scalar unit_cost. A profile carrying a vector unit_cost is rejected, because this function orders costs from the lowest to the highest stratum while n_alloc() orders them by frame row, so a vector written for one would be silently misapplied by the other.

Value

A svyplan_strata object with components:

boundaries: Numeric vector of cutpoints (length n_strata - 1).
n_strata: Number of strata.
n: Total sample size.
cv: Coefficient of variation achieved by the integer allocation in $strata$n (the continuous optimum's cv is kept in params$cv_continuous, and the requested target, if any, in params$cv_target).
strata: Data frame with per-stratum summaries: stratum, lower, upper, N, share, sd, mean, n, and take_all. The N, sd, and mean columns are exactly what n_alloc() expects, so the table can be passed straight to it.
method: Algorithm used.
alloc: Allocation method name (character).
params: List of additional parameters.
converged: Logical (for iterative methods). For "kozak", whether the best run stopped on its no-change rule rather than at max_iter, or the boundary sets were enumerated.

Details

Choosing a method

The four methods differ in approach and use case:

cumrootf: Dalenius-Hodges (1959) cumulative root frequency rule. Non-iterative, does not require n or cv. Best for a quick first look at reasonable boundary positions.
geo: Gunning-Horgan (2004) geometric progression. Non-iterative, requires x > 0. Works well for right-skewed positive variables (e.g. income, revenue) where a log-scale spacing is natural.
lh (default): coordinate-wise optimization inspired by Lavallée-Hidiroglou (1988). It starts from quantile boundaries and repeatedly performs bounded one-dimensional minimizations. Requires n or cv. This is a local heuristic, not the published LH recurrence and not a globally optimal algorithm.
kozak: Kozak's (2004) random search, as described by Baillargeon and Rivest (2011, section 3.3). Each iteration moves one randomly chosen boundary by a random number of positions among the distinct values of x, and keeps the move only when it lowers the objective on a feasible partition. A run stops once the boundaries have not changed for a set number of consecutive iterations, and is then reported converged. That rule does not certify a local optimum, and the method returns the best partition found across its runs, with no global-optimality guarantee. Requires n or cv.

The tuning defaults follow the current stratification package, not the 2011 paper, which uses a maximum step of 3, a stopping count of 100, three runs, and enumeration below 1,000 boundary sets. Here the maximum step is a tenth of the distinct values, rounded up and capped at 100, and the stopping count is ten times the step, between 50 and 500. n_restart runs start from each of the cumulative root frequency, geometric (when x > 0) and quantile boundaries. When there are at most 10,000 possible boundary sets they are all evaluated instead, which gives the exact optimum.

In summary, "lh" is the faster iterative heuristic. "kozak" spends more computation on random moves from several starting points, and is exact on small problems. Use "cumrootf" or "geo" when you do not yet have n or cv.

Allocation is controlled by the alloc parameter. Four methods are available:

proportional: $n_h \propto N_h$.
neyman: $n_h \propto N_h S_h$. Minimizes the national CV when unit costs are equal.
optimal: $n_h \propto N_h S_h / \sqrt{c_h}$. Accounts for differential unit costs.
power: Bankier (1988) allocation, $n_h \propto (S_h/|\mu_h|)(N_h|\mu_h|)^q$. The population CV and default allocation measure N * abs(mean) are recomputed for each candidate partition. Stratum means must be nonzero, and x may not take both signs. At alloc_q = 1 this gives Neyman allocation. At alloc_q = 0 it allocates in proportion to population CVs. The rule does not generally equalize stratum CVs. For a custom measure, first construct the strata, then supply an alloc_measure column to n_alloc().

Stratum allocations are rounded to integers using the ORIC method (Cont and Heidari, 2015), which preserves sum(n) = n while minimizing rounding distortion.

What the boundaries are optimal for

Every quantity the search reads comes from x. The stratum standard deviations $S_h$ are the standard deviations of x inside each candidate stratum, and the variance being minimized is that of the estimated total or mean of x. The boundaries returned are therefore optimal for estimating x itself.

Surveys rarely measure x. It is an auxiliary variable already on the frame, and the study variable y is what will be collected. Nothing here models $E(y \mid x)$ or $\mathrm{Var}(y \mid x)$, so for any y other than x these boundaries are a proxy, good in proportion to how closely y tracks x. A frame with household expenditure stratified for a poverty rate is the usual case, and it is a reasonable one. A frame stratified on establishment size for a variable unrelated to size is not.

Two consequences worth planning around. The reported $cv is the CV for x, not for the survey's own indicators, so it bounds what the design achieves only to the extent that y and x share a stratum structure. And with several study variables no single x is optimal for all of them: choose the x closest to the indicator that matters most, or compare the boundary sets a few candidate x produce before committing.

Design effect, response rate, and the boundaries

The variance evaluated here is the one the rest of the package uses, $\mathrm{deff}\sum_h W_h^2S_h^2(1/(n_h r)-1/N_h)$ with response rate $r$, so cv = 0.05 means the same thing to strata_bound() and n_alloc(), and $n is a fielded sample in both. Handing $strata to n_alloc() with the same cv, deff, and resp_rate reproduces this function's continuous total.

A scalar deff scales the variance of every candidate boundary set equally, so it does not move the boundaries. It changes the n a cv target needs and the cv a given n achieves. Expect the cutpoints to shift only by the local search's own tolerance. Boundaries would respond to a design effect that varied across strata, which a scalar argument cannot express.

The reported $cv is that of the integer allocation in $strata$n. Because strata_bound() rounds up in cv mode and n_alloc() reports a continuous optimum, the two totals differ by the rounding, not by the model.

References

Dalenius, T. and Hodges, J. L. (1959). Minimum variance stratification. Journal of the American Statistical Association, 54(285), 88–101.

Lavallée, P. and Hidiroglou, M. (1988). On the stratification of skewed populations. Survey Methodology, 14(1), 33–43.

Kozak, M. (2004). Optimal stratification using random search method in agricultural surveys. Statistics in Transition, 6(5), 797–806.

Baillargeon, S. and Rivest, L.-P. (2011). The construction of stratified designs in R with the package stratification. Survey Methodology, 37(1), 53–65.

Gunning, P. and Horgan, J. M. (2004). A new algorithm for the construction of stratum boundaries in skewed populations. Survey Methodology, 30(2), 159–166.

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

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

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

Examples

set.seed(867)
x <- rlnorm(500, meanlog = 6, sdlog = 1.5)

# Dalenius-Hodges (non-iterative)
strata_bound(x, n_strata = 4, method = "cumrootf", n = 100)
#> Strata boundaries (Dalenius-Hodges, 4 strata)
#> n = 100, cv = 0.0207, allocation: neyman
#> 
#>  stratum  lower upper   N share     sd    mean  n
#>        1 8.6346   600 299 0.598  150.6   222.9 25
#>        2    600  1800 116 0.232  339.6   994.2 21
#>        3   1800  4600  54 0.108  774.1  2731.8 23
#>        4   4600 30185  31 0.062 6786.6 10142.3 31

# LH (default, iterative)
strata_bound(x, n_strata = 4, n = 100)
#> Strata boundaries (LH-inspired coordinate search, 4 strata, converged)
#> n = 100, cv = 0.0168, allocation: neyman
#> 
#>  stratum  lower  upper   N share     sd   mean  n
#>        1 8.6346 422.82 262 0.524  112.8  183.1 17
#>        2 422.82 1172.6 123 0.246  200.6  725.4 14
#>        3 1172.6 2524.2  59 0.118  381.4 1809.6 13
#>        4 2524.2  30185  56 0.112 6053.2 7141.7 56

# Bankier allocation balances population CVs and stratum totals
strata_bound(x, n_strata = 4, n = 100, alloc = "power", alloc_q = 0.5)
#> Strata boundaries (LH-inspired coordinate search, 4 strata, converged)
#> n = 100, cv = 0.0170, allocation: power (alloc_q = 0.50)
#> 
#>  stratum  lower  upper   N share     sd   mean  n
#>        1 8.6346 441.45 266 0.532  116.1  186.9 21
#>        2 441.45 1232.3 122 0.244  208.0  746.9 13
#>        3 1232.3 2524.2  56 0.112  364.8 1841.5 10
#>        4 2524.2  30185  56 0.112 6053.2 7141.7 56

# With take-all stratum
strata_bound(x, n_strata = 3, n = 80, take_all_above = quantile(x, 0.95))
#> Strata boundaries (LH-inspired coordinate search, 3 strata, converged)
#> n = 80, cv = 0.0403, allocation: neyman
#> 
#>  stratum  lower  upper   N share     sd    mean  n
#>        1 8.6346 1110.2 383 0.766  287.4   352.2 29
#>        2 1110.2 5502.7  92 0.184 1107.0  2450.4 26
#>        3 5502.7  30185  25 0.050 7057.1 11343.1 25

# Under a clustered design and imperfect response, on the same scale
# n_alloc() uses
sb <- strata_bound(x, n_strata = 4, cv = 0.08, deff = 1.8,
                   resp_rate = 0.85)
sb
#> Strata boundaries (LH-inspired coordinate search, 4 strata, converged)
#> n = 46, cv = 0.0765, allocation: neyman
#> 
#>  stratum  lower  upper   N share     sd    mean  n
#>        1 8.6346 441.45 266 0.532  116.1   186.9  5
#>        2 441.45 1384.3 131 0.262  248.4   786.4  5
#>        3 1384.3 5020.2  75 0.150  948.7  2555.1 10
#>        4 5020.2  30185  28 0.056 6919.0 10697.5 26

# The strata table hands straight to n_alloc()
n_alloc(sb$strata, cv = 0.08, deff = 1.8, resp_rate = 0.85)
#> Stratum allocation (neyman, 4 strata)
#> field design: n = 46, expected respondents = 39.1, cv = 0.0765, cost = 46
#> continuous optimum: n = 43.82296, cv = 0.0800, se = 103.0237
#> (resp_rate = 0.85, deff = 1.80)
#> design df = 42

# Or carry the design assumptions in a profile
plan <- svyplan(deff = 1.8, resp_rate = 0.85)
strata_bound(x, n_strata = 4, cv = 0.08, plan = plan)
#> Strata boundaries (LH-inspired coordinate search, 4 strata, converged)
#> n = 46, cv = 0.0765, allocation: neyman
#> 
#>  stratum  lower  upper   N share     sd    mean  n
#>        1 8.6346 441.45 266 0.532  116.1   186.9  5
#>        2 441.45 1384.3 131 0.262  248.4   786.4  5
#>        3 1384.3 5020.2  75 0.150  948.7  2555.1 10
#>        4 5020.2  30185  28 0.056 6919.0 10697.5 26