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 random-restart local search inspired by Kozak.
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_aboveis specified). Must be >= 2.- ...
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
norcv. Required for methods"lh"and"kozak".- cv
Target coefficient of variation (relative standard error). For example,
cv = 0.05means the standard error of the estimated population total or mean should be at most 5 percent of the estimate. Specify at most one ofnorcv. Required for methods"lh"and"kozak".- method
Stratification method:
"cumrootf"(Dalenius-Hodges),"geo"(geometric),"lh"(LH-inspired coordinate optimization), or"kozak"(Kozak-inspired random-restart local search). The short method names are retained for compatibility and do not claim exact implementations of the published algorithms. Default"lh".- alloc
Allocation rule:
"proportional","neyman","optimal", or"power"(Bankier compromise). Default"neyman". See Details.- alloc_q
Bankier power parameter, used only when
alloc = "power". Numeric scalar in \([0, 1]\). Atalloc_q = 1the allocation equals Neyman. Atalloc_q = 0it yields near-equal subnational 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. DefaultNULL(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_aboveform 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 thatcvandnmean 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,
nis the sample fielded andn * resp_ratethe number expected to respond, matchingn_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. DefaultNULL(Freedman-Diaconis rule).- max_iter
Positive whole number of maximum iterations. Applies to
method = "lh"and"kozak"only. Supplying it with another method is an error rather than a silent no-op. DefaultNULL(= 200).- n_restart
Positive whole number of random restarts. Applies to
method = "kozak"only. Supplying it with another method is an error rather than a silent no-op. DefaultNULL(= 10 *n_strata).- plan
Optional
svyplan()object providing design defaults. It can supplyalloc,alloc_q,deff,resp_rate, and a scalarunit_cost. A profile carrying a vectorunit_costis rejected, because this function orders costs from the lowest to the highest stratum whilen_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 inparams$cv_continuous, and the requested target, if any, inparams$cv_target).- strata
Data frame with per-stratum summaries:
stratum,lower,upper,N,share,sd,mean,n, andtake_all. TheN,sd, andmeancolumns are exactly whatn_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).
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
norcv. 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
norcv. This is a local heuristic, not the published LH recurrence and not a globally optimal algorithm.kozak: random-restart adjacent-boundary local search inspired by Kozak (2004). Requires
norcv. It can explore more starting points than"lh", but it can still finish at a local optimum and provides no global-optimality guarantee.
In summary, "lh" is the faster iterative heuristic. "kozak" spends more
computation on random restarts and may find a better local solution. 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) compromise, \(n_h \propto S_h N_h^{q}\). The parameter
alloc_qcontrols the trade-off between national precision (alloc_q = 1, equivalent to Neyman) and near-equal subnational CVs (alloc_q = 0).
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.
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.
See also
predict.svyplan_strata to assign new data to strata,
n_alloc() to distribute a sample across an existing set of strata,
which accepts $strata directly, and svyplan() for reusable design
defaults.
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 428.41 263 0.526 113.6 184.0 17
#> 2 428.41 1221.9 124 0.248 207.0 735.5 15
#> 3 1221.9 2641.1 58 0.116 382.0 1844.6 13
#> 4 2641.1 30185 55 0.110 6077.6 7223.7 55
# Bankier power allocation (compromise between national and subnational CVs)
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.0176, allocation: power (alloc_q = 0.50)
#>
#> stratum lower upper N share sd mean n
#> 1 8.6346 396.92 255 0.510 107.5 176.8 11
#> 2 396.92 1176.7 130 0.260 207.4 708.6 15
#> 3 1176.7 2661.2 60 0.120 392.9 1823.3 19
#> 4 2661.2 30185 55 0.110 6077.6 7223.7 55
# 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 1070.6 378 0.756 276.0 342.3 27
#> 2 1070.6 5502.7 97 0.194 1118.9 2380.8 28
#> 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.0779, allocation: neyman
#>
#> stratum lower upper N share sd mean n
#> 1 8.6346 840.54 350 0.700 220.9 292.8 13
#> 2 840.54 3539.1 109 0.218 704.0 1712.1 13
#> 3 3539.1 13717 36 0.072 2388.4 6498.6 15
#> 4 13717 30185 5 0.010 4835.3 24171.6 5
# 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.0779, cost = 46
#> continuous optimum: n = 43.97615, 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.0779, allocation: neyman
#>
#> stratum lower upper N share sd mean n
#> 1 8.6346 840.54 350 0.700 220.9 292.8 13
#> 2 840.54 3539.1 109 0.218 704.0 1712.1 13
#> 3 3539.1 13717 36 0.072 2388.4 6498.6 15
#> 4 13717 30185 5 0.010 4835.3 24171.6 5