joint_frame <- data.frame(
stratum = c("NU", "NR", "SU", "SR"),
region = c("North", "North", "South", "South"),
residence = c("Urban", "Rural", "Urban", "Rural"),
N = c(1000, 2000, 1500, 1000),
unit_cost = c(1, 1.2, 1.5, 1)
)
joint_measures <- data.frame(
stratum = rep(joint_frame$stratum, 2),
name = rep(c("vaccination", "income"), each = 4),
p = c(0.50, 0.40, 0.60, 0.30, rep(NA, 4)),
mean = c(rep(NA, 4), 50, 55, 60, 45),
sd = c(rep(NA, 4), 10, 12, 15, 9)
)
joint_targets <- data.frame(
name = c("vaccination", "vaccination", "income"),
domain = c(".overall", "region", "residence"),
level = c(NA, "North", "Urban"),
cv = c(0.05, 0.08, NA),
moe = c(NA, NA, 2)
)4 Joint allocation across indicators and domains
A joint allocation assigns one sample across strata while satisfying several indicator and domain requirements. Use this interface when separate single-indicator allocations would compete for the same fieldwork sample.
Chapter 3 introduces strata, costs, and operational rounding for a single indicator. Chapter 2 shows how several requirements determine a total before allocation.
4.1 Three input tables
The joint interface separates population structure, anticipated measurements, and publication requirements:
| Table | One row represents | Main contents |
|---|---|---|
frame |
An atomic stratum | Population size, domain labels, costs, and design bounds. |
measures |
One indicator in one stratum | Anticipated proportion or mean and standard deviation. |
targets |
One required estimate | Indicator, domain and level, and precision target. |
For element-level designs, a proportion in measures$p uses the finite-population variance N * p * (1 - p) / (N - 1) for its stratum, with N > 1. This matches the convention in prec_prop(). Explicit variances and standard deviations are used as supplied. The multistage designs introduced below retain the working proportion variance p * (1 - p).
An atomic stratum has one value for every domain classification. Crossing region with residence therefore produces four frame rows:
The targets request an overall vaccination CV, a separate CV for vaccination in the North, and a margin of error for urban income. The optimizer minimizes variable cost while meeting all three:
joint_fit <- n_alloc(
joint_frame,
measures = joint_measures,
targets = joint_targets,
min_n_stratum = 2
)
joint_fit
#> Joint constrained allocation (Bethel)
#> question: cheapest design meeting every precision target
#> field design: n = 425, cost = 508 (3 targets, all pass)
#> continuous optimum: n = 425.113, cost = 508 (integerizing costs +0.02%)
#> binding: vaccination@.overall:cv (target 0.05, achieved 0.05)
joint_fit$constraints
#> constraint name domain level
#> 1 vaccination@.overall:cv vaccination .overall <NA>
#> 2 vaccination@region=North:cv vaccination region North
#> 3 income@residence=Urban:moe income residence Urban
#> .metric .target .achieved .ratio .residual .tolerance
#> 1 cv 0.05 0.05000000 1.0000000 2.459284e-10 1e-06
#> 2 cv 0.08 0.07030389 0.8787986 -1.212014e-01 1e-06
#> 3 moe 2.00 1.85266616 0.9263331 -7.366692e-02 1e-06
#> .se .cv .moe .rmoe .pass .binding
#> 1 0.02272727 0.05000000 0.04454464 0.09799820 TRUE TRUE
#> 2 0.03046502 0.07030389 0.05971034 0.13779309 TRUE FALSE
#> 3 0.94525521 0.01687956 1.85266616 0.03308332 TRUE FALSE
#> .multiplier .sensitivity
#> 1 0.03001964 -18762.27
#> 2 0.00000000 0.00
#> 3 0.00000000 0.00$constraints shows the target and achieved precision for every requirement. Requirements closest to their limits determine the allocation. The remaining requirements can be better than requested because they share the same sample.
4.2 Continuous and fieldwork allocations
The continuous optimum and the whole-unit recommendation are both retained:
joint_fit$detail[, c("stratum", "n", "n_int", "weight")]
#> stratum n n_int weight
#> 1 NU 86.67423 87 11.53746
#> 2 NR 155.00879 156 12.90249
#> 3 SU 103.99173 104 14.42422
#> 4 SR 79.43824 78 12.58840The integer allocation is constructed to remain feasible. Check both designs directly:
prec_alloc(joint_fit)
#> Joint allocation precision (3 constraints)
#> targets: all pass
#> showing 1 binding of 3
#>
#> constraint .metric .target .achieved .pass
#> vaccination@.overall:cv cv 0.05 0.05 TRUE
#>
#> ... see $detail for all rows
prec_alloc(joint_fit, n = joint_fit$detail$n_int)
#> Joint allocation precision (3 constraints)
#> targets: all passA modified named allocation can be passed in any row order because names are matched to frame$stratum. When comparing a simulation with the plan, compare the empirical precision with $constraints$.achieved, not only with the requested target, because integerization can make the fieldwork design slightly more precise.
4.3 Fixed-take multistage allocation
The same interface supports two-stage designs when the number of ultimate units sampled per primary sampling unit (PSU) is fixed. The frame holds PSU populations, takes, and stage costs. The measures table holds the indicator-specific PSU homogeneity. Only the number of PSUs is optimized.
cluster_frame <- within(joint_frame, {
unit_cost <- NULL
N_psu <- c(100, 160, 120, 90)
n_per_psu <- c(8, 10, 12, 7)
cost_psu <- c(300, 400, 450, 350)
cost_ssu <- c(25, 30, 35, 28)
})
cluster_measures <- transform(
joint_measures,
icc_psu = rep(c(0.03, 0.05, 0.08, 0.04), 2),
var_ratio_psu = 1
)cluster_fit <- n_alloc(
cluster_frame,
measures = cluster_measures,
targets = joint_targets,
min_n_stratum = 20
)
cluster_fit$detail[, c(
"stratum", "n_psu_int", "n_per_psu", "n_int", "weight"
)]
#> stratum n_psu_int n_per_psu n_int weight
#> 1 NU 14 8 112 9.244275
#> 2 NR 22 10 220 9.121262
#> 3 SU 16 12 192 8.152297
#> 4 SR 12 7 84 11.130697
prec_alloc(cluster_fit, n = cluster_fit$detail$n_int)
#> Joint allocation precision (3 constraints)
#> targets: all passPublic n remains the ultimate-unit count. The operational identity is n_int = n_psu_int * n_per_psu. A three-stage frame adds N_ssu, n_per_ssu, cost_tsu, and the corresponding indicator-specific SSU homogeneity. take_all is unavailable in fixed-take multistage mode because a census of PSUs does not imply a census of ultimate units.
4.4 Certainty PSUs from a register
For a two-stage design, a PSU register can replace N_psu. It has one row per PSU with its stratum and size, and its stratum totals must match frame$N. Every PSU must hold at least its stratum’s take, so merge smaller ones with a neighbour before planning. n_alloc() identifies PSUs whose implied probability-proportional-to-size inclusion probability reaches 1. That includes the whole number of PSUs the remainder is drawn as: rounding it up can push a PSU below the size threshold to probability 1, and $psu$.certainty_source then reads "operational". A certainty column can mark additional certainty PSUs.
certainty_frame <- cluster_frame
certainty_frame$N_psu <- NULL
psu_register <- do.call(rbind, lapply(seq_len(nrow(certainty_frame)), function(i) {
population <- certainty_frame$N[i]
size <- rep(12, 20)
size[1] <- floor(population / 3)
size[2] <- population - sum(size[-2])
data.frame(stratum = certainty_frame$stratum[i], N = size)
}))The remainder’s PSUs are selected with probability proportional to size, so its icc_psu is the size-weighted share of the stratum variance that lies between PSUs. varcomp() returns it when prob is each PSU’s share of its stratum’s units. Use every PSU of the stratum, since which ones end up certainty depends on the allocation. Here census stands for a census or pilot with one row per unit:
set.seed(2)
psu_register$psu_id <- seq_len(nrow(psu_register))
census <- data.frame(
stratum = rep(psu_register$stratum, psu_register$N),
psu_id = rep(psu_register$psu_id, psu_register$N)
)
psu_effect <- rnorm(nrow(psu_register))
census$vaccination <- rbinom(
nrow(census), 1, plogis(-0.2 + 0.5 * psu_effect[census$psu_id])
)
census$income <- 50 + 3 * psu_effect[census$psu_id] +
rnorm(nrow(census), 0, 11)
census$share <- ave(census$psu_id, census$psu_id, FUN = length) /
ave(census$psu_id, census$stratum, FUN = length)
icc_of <- function(name) {
vc <- varcomp(reformulate("psu_id", name), strata = ~stratum,
prob = ~share, data = census)
data.frame(stratum = vc$strata$stratum, name = name,
icc_psu = vc$strata$icc_psu)
}
certainty_measures <- merge(
joint_measures, rbind(icc_of("vaccination"), icc_of("income"))
)
certainty_measures
#> stratum name p mean sd icc_psu
#> 1 NR income NA 55 12 0.16790848
#> 2 NR vaccination 0.4 NA NA 0.13327002
#> 3 NU income NA 50 10 0.05051815
#> 4 NU vaccination 0.5 NA NA 0.04113408
#> 5 SR income NA 45 9 0.20123765
#> 6 SR vaccination 0.3 NA NA 0.14481297
#> 7 SU income NA 60 15 0.07217228
#> 8 SU vaccination 0.6 NA NA 0.05855218Without prob, varcomp() assumes PSUs of equal chance and measures the spread of PSU totals, which carries the spread of PSU sizes. On a register as unequal as this one that is not the homogeneity the model needs:
varcomp(income ~ psu_id, strata = ~stratum, data = census)$strata$icc_psu
#> [1] 0.9467741 0.9525848 0.9335702 0.9546080varcomp() needs at least two PSUs in a stratum. A stratum with a single PSU has no variance between PSUs, so leave its icc_psu as NA and n_alloc() uses 0.
certainty_fit <- n_alloc(
certainty_frame,
measures = certainty_measures,
targets = joint_targets,
psu = psu_register,
min_n_stratum = 20
)
certainty_fit$detail[, c(
"stratum", "n_psu_certain", "n_psu_draw", "n_int"
)]
#> stratum n_psu_certain n_psu_draw n_int
#> 1 NU 2 3 95
#> 2 NR 2 2 180
#> 3 SU 2 2 129
#> 4 SR 2 3 85n_psu_certain counts the certainty PSUs and n_psu_draw the remainder PSUs the operational design draws, so n_int = n_certain_int + n_psu_draw * n_per_psu. $psu records each PSU’s classification, threshold, and source. The certainty takes are rounded up one PSU at a time, so $operational$constraints assesses the field design on its own takes, and the remainder gains whole PSUs where that falls short of a target. prec_alloc(certainty_fit) reproduces the continuous plan’s precision. Under stage costs a certainty PSU’s visit is a fixed cost, and a budget is refused with a register.
Practitioners often relax the rule and take with certainty any PSU whose inclusion probability reaches, say, 0.80. certainty_cutoff sets that level, for every stratum or per stratum through a frame column, and predict(certainty_fit, data.frame(certainty_cutoff = c(1, 0.9, 0.8))) prices the choice. Here it changes nothing until the cutoff falls to 0.1, because no PSU sits between the two dominant ones and the small ones. The register example under ?n_alloc shows a register where it does.
Area samples often go further and cut each stratum’s remainder into zones of about equal size, drawing two PSUs from each so that every stratum’s variance can be estimated. n_psu_per_zone = 2 plans that design, and an optional zone_order column in the register sets the order the zones are cut in. n_psu_per_zone = 1 draws one PSU per zone, and $psu$.pair groups the zones in advance for a collapsed variance estimate. A zone holding a single PSU is a census, which can lower the sample, though not on this register, whose small PSUs are all of one size. predict() compares the settings:
predict(certainty_fit, data.frame(n_psu_per_zone = c(NA, 1, 2)))
#> n_psu_per_zone n_psu_certain n_psu_draw n_zone n
#> 1 NA 8 10 0 467.8852
#> 2 1 8 10 10 467.8852
#> 3 2 8 12 6 467.8852
#> cost n_int cost_int .feasible
#> 1 20041.55 489 21320 TRUE
#> 2 20041.55 489 21320 TRUE
#> 3 20041.55 504 22366 TRUESome planning tools also move small groups of remainder PSUs into the certainty part, so that each group yields at least two PSUs for variance estimation. That grouping is a selection decision, and n_alloc() does not make it. To plan with it, form the groups from the PSU sizes and $psu$.threshold, set certainty = TRUE in the register for the PSUs to promote, and solve again. Those PSUs report .certainty_source as "supplied", and the printed plan counts them. Expect the certainty count to differ from such a tool’s until you do.
4.5 Explore controls that remain movable
For a fixed-take multistage result, predict() can re-solve the stored problem over candidate takes:
predict(cluster_fit, data.frame(n_per_psu = c(6, 10, 14, 20)))
#> n_per_psu n_psu n cost n_psu_int n_int cost_int
#> 1 6 87.43080 524.5848 49366.72 88 528 49752
#> 2 10 60.14504 601.4504 41283.14 60 600 41350
#> 3 14 48.25956 675.6338 39003.24 49 686 39614
#> 4 20 39.11771 782.3542 38764.28 40 800 39750
#> .feasible
#> 1 TRUE
#> 2 TRUE
#> 3 TRUE
#> 4 TRUEChanging precision targets or the measures table defines a different planning problem. Edit those inputs and call n_alloc() again. See ?predict.svyplan for the supported controls and output columns for each joint-allocation mode.
Chapter 5 plans a single indicator with independently chosen cluster counts and takes.