---
title: "An area frame with points, a list frame, and their joint design"
output:
  rmarkdown::html_vignette:
    css: fieldopt.css
vignette: >
  %\VignetteIndexEntry{An area frame with points, a list frame, and their joint design}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- Generated from area-and-dual-frame.Rmd.orig by tools/precompute-vignettes.R; edit the .orig file. -->



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


``` r
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.


``` r
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.


``` r
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.


``` r
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.


``` r
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:


``` r
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.


``` r
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.


``` r
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.


``` r
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.


``` r
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.


``` r
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.


``` r
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).


``` r
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.