An area frame with points, a list frame, and their joint design

In fieldopt: Area-Frame and Dual-Frame Survey Design under Real Field Cost

Agricultural 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        236

Each 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.9

Stratified 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] 675

In 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            164

The 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.9

The 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          3

What 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     3

The 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.0

Dual-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.7470

Auxiliaries: 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.275

References

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.

fieldopt documentation built on Oct. 22, 2026, 1:07 a.m.