An area frame with points, a list frame, and their joint design
Source:vignettes/area-and-dual-frame.Rmd
area-and-dual-frame.RmdAgricultural surveys often combine an area frame (a
grid of square cells classified by intensity of agricultural use from
satellite imagery, with sample points inside the selected cells
identifying the fields and so the establishments) and a list
frame (large and specialised establishments, cheap to reach).
This vignette walks through such a design with fieldopt:
stratified selection of cells, points inside them, the multiplicity
estimator, the cost of the field work, and the dual-frame and two-stage
allocations that put the pieces together. The data are synthetic; the
structure is the usual one of area frames built on a regular grid
(square segments of about 1 km²) with a list of large producers
alongside.
The area frame: cells with intensity classes
library(fieldopt)
library(dplyr, warn.conflicts = FALSE)
library(tidyr)
set.seed(2026)
side <- 30
cells <- expand_grid(x = seq_len(side), y = seq_len(side)) |>
mutate(
unit = paste0("c", row_number()),
use = 0.5 * x / side + 0.3 * sin(y / 4) + rnorm(n(), sd = 0.15),
intensity = cut(use, c(-Inf, 0.25, 0.5, Inf), c("low", "mid", "high")),
size = c(low = 1, mid = 2, high = 4)[as.character(intensity)]
) |>
select(-use)
cells |> count(intensity)
#> # A tibble: 3 × 2
#> intensity n
#> <fct> <int>
#> 1 low 380
#> 2 mid 284
#> 3 high 236Each cell holds a few establishments; the agricultural area of an establishment inside a cell is what the point sampling needs. Here every cell has two establishments with 40 % and 20 % of its area, and a crop value per establishment that follows the intensity.
areas <- bind_rows(
cells |> transmute(unit, establishment = paste0(unit, "-1"), area = 0.4),
cells |> transmute(unit, establishment = paste0(unit, "-2"), area = 0.2)
) |>
left_join(cells |> select(unit, size), by = "unit") |>
mutate(y = round(rgamma(n(), shape = 4, rate = 4 / (20 * area * size)), 1)) |>
select(-size)
truth <- sum(areas$y)
truth
#> [1] 22447.9Stratified selection of cells and points inside them
Cells are selected within intensity strata, more where use is
intense, with the local pivotal method so that the sample spreads over
the territory; then 3, 6 or 9 points per cell depending on the stratum,
laid out as a regular grid with one random offset per cell
(layout = "systematic"), which keeps each point uniform
over the cell while spreading them.
n_h <- c(low = 20, mid = 35, high = 45)
cells_s <- select_units(cells, n = n_h, strata = "intensity", seed = 7)
cells_s
#> Sample of 100 units from 900 by "lpm", 3 strata.
#> # A tibble: 900 × 7
#> x y unit intensity size pi sampled
#> * <int> <int> <chr> <fct> <dbl> <dbl> <lgl>
#> 1 1 1 c1 low 1 0.0526 FALSE
#> 2 1 2 c2 low 1 0.0526 FALSE
#> 3 1 3 c3 low 1 0.0526 FALSE
#> 4 1 4 c4 mid 2 0.123 FALSE
#> 5 1 5 c5 low 1 0.0526 FALSE
#> 6 1 6 c6 low 1 0.0526 FALSE
#> 7 1 7 c7 low 1 0.0526 FALSE
#> 8 1 8 c8 low 1 0.0526 FALSE
#> 9 1 9 c9 mid 2 0.123 TRUE
#> 10 1 10 c10 low 1 0.0526 FALSE
#> # ℹ 890 more rows
ppc <- c(low = 3, mid = 6, high = 9)
points <- cells_s |>
select_points(points_per_cell = ppc, cell_size = 1, layout = "systematic", seed = 7)
nrow(points)
#> [1] 675In the field each point is visited and the establishment it falls on is interviewed about its whole operation. Here the field is simulated: a point hits establishment 1 with probability 0.4, establishment 2 with probability 0.2, and otherwise falls on non-agricultural land.
set.seed(8)
hits <- points |>
mutate(u = runif(n()), which = case_when(u < 0.4 ~ "-1", u < 0.6 ~ "-2")) |>
filter(!is.na(which)) |>
transmute(unit, establishment = paste0(unit, which)) |>
left_join(areas |> select(establishment, y), by = "establishment")
c(hits = nrow(hits), establishments = n_distinct(hits$establishment))
#> hits establishments
#> 401 164The multiplicity estimator
An establishment can be hit several times; the estimator divides the value of every hit by the expected number of hits of that establishment under the design, which depends on the cells it spans, their inclusion probabilities and their point densities.
eh <- expected_hits(areas, cells_s, points_per_cell = ppc, cell_size = 1)
hits <- hits |> left_join(eh |> select(establishment, expected_hits), by = "establishment")
point_estimator(hits, cells_s)
#> # A tibble: 1 × 11
#> total variance se cv variance_srs n N variance_method df
#> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int> <chr> <int>
#> 1 22034. 1178527. 1086. 0.0493 3243984. 100 900 stratified local-… 97
#> # ℹ 2 more variables: n_hits <int>, n_establishments <int>
truth
#> [1] 22447.9The variance uses the cell design (local-mean estimator on the cell contributions). The same cells selected as a spatially ordered systematic sample with four interpenetrating replicates give a replicate-based variance instead:
cells <- cells |> left_join(areas |> group_by(unit) |> summarise(crop = sum(y)), by = "unit")
cells_r <- cells |>
select_units(
n = n_h, strata = "intensity", method = "systematic", replicates = 4, seed = 7
)
design_variance(cells_r, "crop")
#> # A tibble: 1 × 9
#> total variance se cv variance_srs n N variance_method df
#> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int> <chr> <int>
#> 1 22562. 517031. 719. 0.0319 2488872. 92 900 replicates 3What the field work costs
The cells are visited from a depot; the routing gives the travel, the cost model the money. The cost per visited cell is the input the allocations below need, and it depends on how many cells are visited and how spread they are, which is why it is computed rather than assumed.
frame <- bind_rows(
tibble(unit = "depot", x = side / 2, y = side / 2),
cells |> select(unit, x, y)
)
model <- field_cost_model(
per_travel = 2.5, per_unit = 60, per_interview = 25,
interviews_per_unit = 4, currency = "BRL"
)
routed_unit_cost(
frame, "depot", model,
n = sum(n_h), method = "euclidean",
n_rep = 3, iterations = 60
)
#> # A tibble: 1 × 6
#> n cost_per_unit travel_per_unit cost_mean routes_mean n_rep
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 100 167. 2.67 16667. 1 3The design effect of point sampling
The dual-frame allocation below needs the design effect of the area
sample relative to simple random sampling of establishments.
point_design_effect() simulates the whole area design
(cells, points, hits, multiplicity estimator) and returns the bias, the
empirical variance, the mean of the variance estimates, the number of
interviews and the design effect.
deff <- point_design_effect(
cells, areas,
n = n_h, strata = "intensity",
points_per_cell = ppc, cell_size = 1, n_sim = 60, seed = 9
)
deff
#> # A tibble: 1 × 10
#> n_sim total mean_estimate bias variance cv mean_variance_estimate
#> * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 60 22448. 22355. -93.1 1554792. 0.0555 1349171.
#> # ℹ 3 more variables: interviews <dbl>, variance_srs <dbl>, deff <dbl>The bias is negligible and the variance estimator of
point_estimator() tracks the empirical variance; the design
effect is what enters deff_a.
Two-stage allocation: how many cells, how many points
With the between-cell and within-cell variance components (from a
pilot or from the frame) and the costs per cell and per point, Cochran’s
rule gives the number of points per cell and then the number of cells
for a target precision. The cost per cell is not a constant: a route
shares its travel among the cells it visits, so the routed cost of
n cells grows like c0 + a n + b sqrt(n).
two_stage_design() routes samples at a few sizes, fits that
curve and allocates with the marginal cost of a cell,
which is what Cochran’s rule needs when the cost is not linear (the
average cost per cell would give too few cells), iterating until the
number of cells settles. This is the travel-cost model of Hansen,
Hurwitz and Madow (1953), whose term C0 sqrt(m) they
derived from points on a regular grid; here b plays the
role of C0 and is measured by routing the frame itself. The
history shows the marginal cost used in each round next to the average,
and routed_cost_curve() gives the curve on its own.
ts <- two_stage_design(
frame, "depot", model,
m_secondary = 9,
s2_between = var(cells$crop / 9), s2_within = var(areas$y) * 2,
target_cv = 0.05, mean = mean(cells$crop) / 9,
method = "euclidean", n_rep = 3, iterations = 60
)
ts
#>
#> ── Two-stage allocation ────────────────────────────────────────────────────────
#> n = 158 primary units with m = 9 secondary units each (optimal m 9): cost
#> 45860, variance of the total 1256000 (SE 1121, CV 4.99%).
ts$history
#> # A tibble: 2 × 7
#> round c1 c1_average n m cost variance_total
#> <int> <dbl> <dbl> <int> <int> <dbl> <dbl>
#> 1 1 90.1 65.3 158 9 45863. 1255934.
#> 2 2 62.4 65.3 158 9 45863. 1255934.
ts$cost_curve
#> # A tibble: 5 × 3
#> n cost_mean cost_sd
#> <int> <dbl> <dbl>
#> 1 95 6360. 4.19
#> 2 126 8310. 12.0
#> 3 158 10312. 11.8
#> 4 198 12805. 11.7
#> 5 253 16209. 12.0Dual-frame allocation: area plus list
The list frame holds the large establishments, which are also on the
area frame (the overlap ab); small establishments are on
the area frame only (a); a few specialised units are on the
list only (b). Given the sizes, means and standard
deviations of the three domains, the cost per list unit and the design
effect from the simulation, dual_frame_design() runs
Hartley’s allocation with the marginal routed cost of an area unit, in
the same way, iterating until the area sample size settles;
theta = "screening" gives the design in which the area
sample does not use the overlap.
domains <- tibble(
domain = c("a", "ab", "b"),
size = c(1500, 250, 50),
mean = c(6, 30, 45), sd = c(5, 20, 35)
)
df <- dual_frame_design(
frame, "depot", model, domains,
cost_b = 45, deff_a = deff$deff,
target_cv = 0.05, method = "euclidean", n_rep = 3, iterations = 60
)
df
#>
#> ── Dual-frame allocation ───────────────────────────────────────────────────────
#> Minimum cost for target variance 878900: n_A = 75 (frame A, cost 163.2/unit),
#> n_B = 129 (frame B, cost 45/unit), theta = 0.163 (optimised).
#> Cost 18390; variance 874000 (frame A 679000, frame B 195000); CV 4.99% of the
#> total 18750.
#> Expected overlap units: 10.7 in the A sample, 107.5 in the B sample.
df$history
#> # A tibble: 3 × 8
#> round cost_a cost_a_average n_a n_b theta cost variance
#> <int> <dbl> <dbl> <int> <int> <dbl> <dbl> <dbl>
#> 1 1 190. 168. 73 136 0.161 18375. 875540.
#> 2 2 163. 168. 75 129 0.163 18386. 874037.
#> 3 3 163. 168. 75 129 0.163 18386. 873980.
tidy(df)
#> # A tibble: 2 × 8
#> frame n cost_per_unit cost deff variance overlap_units
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 A 75 163. 12239. 0.753 679399. 10.7
#> 2 B 129 45 5805 1 194581. 108.
#> # ℹ 1 more variable: weight_on_overlap <dbl>
dual_frame_allocation(
domains,
cost_a = df$cost_a, cost_b = 45, deff_a = deff$deff,
theta = "screening", target_cv = 0.05
)
#>
#> ── Dual-frame allocation ───────────────────────────────────────────────────────
#> Minimum cost for target variance 878900: n_A = 82 (frame A, cost 163.2/unit),
#> n_B = 144 (frame B, cost 45/unit), theta = 0 (fixed).
#> Cost 19860; variance 876100 (frame A 692000, frame B 184000); CV 4.99% of the
#> total 18750.
#> Expected overlap units: 11.7 in the A sample, 120 in the B sample.Estimating from both frames
The domains come from the overlap between the frames:
frame_overlap() places each listed establishment in the
cell it falls in (by coordinates, by polygons, or by a table of pairs
from record linkage, for example with the reclin2 package)
and marks it ab or b, and tells each cell
which listed establishments it holds, which is what the teams use to
screen the overlap out of the area sample.
Once both samples are in, dual_frame_estimator()
combines them. Each sampled unit carries its domain (a or
ab on the area frame, b or ab on
the list); Hartley’s estimator mixes the two estimates of the overlap
with a weight theta that can be fixed, set to zero
(screening) or estimated to minimise the variance, and the
Fuller-Burmeister estimator adds the difference between the two
estimates of the overlap size. The variances come from the variance
estimator of each design.
cells <- cells |>
mutate(
domain = if_else(intensity == "high" & runif(n()) < 0.6, "ab", "a"),
crop_est = if_else(domain == "ab", 60, 12) + rnorm(n(), sd = 3)
)
overlap <- cells |> filter(domain == "ab") # the same holdings, seen from the list
list_frame <- bind_rows(
overlap |> transmute(unit = paste0("l", row_number()), domain, crop_est),
tibble(unit = paste0("b", 1:12), domain = "b", crop_est = 90 + rnorm(12, sd = 5))
) |>
mutate(x = runif(n()), y = runif(n()))
sa <- cells |> select_units(n = 40, strata = "intensity", seed = 3)
sb <- list_frame |> select_units(n = 20, method = "srs", seed = 4)
dual_frame_estimator(sa, sb, y = "crop_est", domain = "domain")
#> # A tibble: 1 × 13
#> total variance se cv estimator theta theta_fixed y_a y_ab_a y_ab_b
#> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <lgl> <dbl> <dbl> <dbl>
#> 1 20078. 186923. 432. 0.0215 hartley 0.217 FALSE 9797. 10005. 9470.
#> # ℹ 3 more variables: y_b <dbl>, n_a <int>, n_b <int>
dual_frame_estimator(sa, sb, y = "crop_est", domain = "domain", theta = "screening")
#> # A tibble: 1 × 13
#> total variance se cv estimator theta theta_fixed y_a y_ab_a y_ab_b
#> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <lgl> <dbl> <dbl> <dbl>
#> 1 19962. 425106. 652. 0.0327 hartley 0 TRUE 9797. 10005. 9470.
#> # ℹ 3 more variables: y_b <dbl>, n_a <int>, n_b <int>
dual_frame_estimator(
sa, sb,
y = "crop_est", domain = "domain", estimator = "fuller-burmeister"
)
#> # A tibble: 1 × 13
#> total variance se cv estimator beta_1 beta_2 y_a y_ab_a y_ab_b y_b
#> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 20077. 186851. 432. 0.0215 fuller-b… 0.253 -2.26 9797. 10005. 9470. 695.
#> # ℹ 2 more variables: n_a <int>, n_b <int>The selection scales to national frames: the local pivotal method and the systematic and cube methods run on a million cells in about a second (a grid index keeps the nearest-neighbour searches local), and the selection is done within strata, so the size of the whole frame is not a constraint.
For everything else the estimation needs, calibration to known
totals, domain estimates, nonresponse adjustment,
as_svydesign() hands a sample to the survey
package as a design object with the right weights, strata and, for
replicated systematic samples, the replicate weights.
d <- as_svydesign(sa)
survey::svytotal(~crop_est, d)
#> total SE
#> crop_est 19802 1632.6
survey::svyby(~crop_est, ~intensity, d, survey::svytotal)
#> intensity crop_est se
#> low low 4873.614 212.3148
#> mid mid 3815.881 160.8595
#> high high 11112.423 1610.7470Auxiliaries: the ratio estimator and balanced samples
Area frames come with auxiliaries known for every cell: the area of
farmland from the land-cover map, the number of holdings from the last
census. The ratio estimator uses one such total at the estimation stage,
and the cube method uses several at the selection stage, drawing a
sample whose Horvitz-Thompson estimates of the auxiliaries match their
totals. spatial_balance() measures how evenly a sample
covers the territory (Stevens and Olsen’s Voronoi measure: zero is
perfect).
cells <- cells |>
mutate(
farmland = size * runif(n(), 15, 35),
crop_r = 0.5 * farmland + rnorm(n(), sd = 3)
)
s_lpm <- cells |> select_units(n = 40, strata = "intensity", seed = 5)
ratio_estimator(s_lpm, y = "crop_r", x = "farmland") |>
select(total, se, variance, variance_ht)
#> # A tibble: 1 × 4
#> total se variance variance_ht
#> <dbl> <dbl> <dbl> <dbl>
#> 1 23722. 411. 169014. 956341.
s_cube <- cells |>
select_units(n = 40, strata = "intensity", method = "cube", balance = "farmland", seed = 5)
s_srs <- cells |> select_units(n = 40, strata = "intensity", method = "srs", seed = 5)
tibble(
design = c("lpm", "cube", "srs"),
balance = c(spatial_balance(s_lpm), spatial_balance(s_cube), spatial_balance(s_srs))
)
#> # A tibble: 3 × 2
#> design balance
#> <chr> <dbl>
#> 1 lpm 0.129
#> 2 cube 0.368
#> 3 srs 0.275References
Cochran, W. G. (1977). Sampling Techniques, third edition. Wiley.
Hartley, H. O. (1974). Multiple frame methodology and selected applications. Sankhya C, 36, 99–118.
Lohr, S. L. and Rao, J. N. K. (2006). Estimation in multiple-frame surveys. Journal of the American Statistical Association, 101(475), 1019–1030.
Nealon, J. P. (1984). Review of the multiple and area frame estimators. USDA Statistical Reporting Service, Staff Report 80.
Deville, J.-C. and Tillé, Y. (2004). Efficient balanced sampling: the cube method. Biometrika, 91(4), 893–912.
Fuller, W. A. and Burmeister, L. F. (1972). Estimators for samples selected from two overlapping frames. Proceedings of the Social Statistics Section, American Statistical Association, 245–249.
Stevens, D. L. and Olsen, A. R. (2004). Spatially balanced sampling of natural resources. Journal of the American Statistical Association, 99, 262–278.