Adaptive cluster sampling in R

adaptive sampling

survey design

abundance

ecology tutorial

Build adaptive cluster sampling from scratch in R for rare, clustered species: networks, edge units, why the naive mean is biased, and the Hansen-Hurwitz fix.

Author

Tidy Ecology

Published

2026-07-26

Rare, clustered species are the hardest thing a fixed design has to cope with. Most quadrats hold nothing, a handful hold almost everything, and a simple random sample of a few dozen units mostly counts zeros. Adaptive cluster sampling turns that structure into an advantage: when a selected unit contains the species, you sample its neighbours too, and keep expanding as long as you keep finding it. The sample follows the organism.

The catch is that the mechanism which makes it efficient also breaks the obvious estimator. If you only expand where you find animals, the units you end up observing over-represent the dense patches, so the plain average of everything you sampled runs high. This post builds the design from scratch on a synthetic grid, shows exactly how biased the naive mean is, and derives the Hansen-Hurwitz estimator that puts it right.

A population built to defeat a fixed design

The frame is a grid of 400 quadrats. The species sits in a few tight clusters on an otherwise empty background: 87% of the quadrats are empty, and the occupancy (the fraction with at least one individual) is only 13.0%. The true quadrat mean is 0.295 individuals and the true total is 118. On a real survey you would not know these; here they are the yardstick.

A 20 by 20 grid of quadrats shaded by count. Most cells are pale (zero), with four small dark clusters of high counts scattered across the grid. — Figure 1: The synthetic population. Counts are concentrated in a few clusters on an empty background, the situation adaptive cluster sampling is built for.

Under simple random sampling you would draw a fixed number of quadrats and stop. With 13% occupancy most of them are zeros, so the estimate of the mean is dominated by whether you happen to land on a cluster. Adaptive cluster sampling changes the second step: it keeps sampling around any quadrat where the species turns up.

The design: condition, neighbourhood, networks

Three ingredients define the design.

The condition decides when to expand. Here a quadrat “satisfies the condition” if its count is at least one, that is, the species is present. The neighbourhood decides where to expand: we use the four rook neighbours (up, down, left, right). The rule is then recursive. Draw an initial simple random sample; whenever a selected quadrat satisfies the condition, add its neighbours; whenever a newly added quadrat also satisfies the condition, add its neighbours too; stop when no new occupied quadrats appear.

The units partition into networks. Connected runs of occupied quadrats form one network each: selecting any one of them pulls in all of them. Every quadrat that does not satisfy the condition is a network on its own. The occupied quadrats that sit on the boundary of a network, its non-satisfying neighbours, are the edge units: they get observed because they are adjacent to a network, but they do not trigger further expansion. build_networks() above finds these by a breadth-first search over the occupied cells.

On this population there are 4 occupied networks, with sizes 19, 12, 12, 9 quadrats. Everything else is a singleton.

## one adaptive draw from an initial sample of n1 = 30 quadrats
set.seed(7)
n1 <- 30L
d1 <- acs_draw(pop, nets, ubn, n1)
c(initial = n1, distinct_networks_hit = length(d1$hit),
  final_sample = d1$n_final, naive_mean = round(d1$naive, 3))

              initial distinct_networks_hit          final_sample 
               30.000                28.000                86.000 
           naive_mean 
                0.767

Thirty initial quadrats expanded to 86 observed quadrats, because a couple of them landed on clusters and dragged the whole cluster (plus its edge) into the sample. The plain average of everything observed in that draw was 0.767, well above the truth of 0.295.

The grid with the initial 30 sampled quadrats marked as open squares, scattered at random, and the adaptively added quadrats marked as filled squares, forming compact blocks around two of the clusters. — Figure 2: One adaptive draw. Open squares are the initial simple random sample; filled squares are quadrats added by expansion. Where an initial unit lands on a cluster, the whole network and its edge come in.

Why the naive mean is biased

The bias is not a fluke of one draw. Average the naive mean over many initial samples and it stays high, because expansion is triggered by high values: you always add neighbours around counts, never around empties, so the observed set is enriched for occupied quadrats relative to the frame. Formally the naive mean is not the mean of a probability sample of the frame at all, and it has no reason to be unbiased.

mc <- function(pop, nets, ubn, n1, B, seed) {
  set.seed(seed); naive <- hh <- nfin <- numeric(B); cov_hh <- logical(B)
  mu <- mean(pop$y)
  for (b in 1:B) {
    d <- acs_draw(pop, nets, ubn, n1)
    naive[b] <- d$naive; hh[b] <- d$hh; nfin[b] <- d$n_final
    se <- sqrt(max(hh_var(pop, nets, d$init, n1), 0))
    cov_hh[b] <- abs(d$hh - mu) <= 1.96 * se
  }
  list(naive = naive, hh = hh, nfin = nfin, cov_hh = cov_hh)
}
M <- mc(pop, nets, ubn, n1 = 30L, B = 20000L, seed = 202L)
c(true_mu = round(mu, 4), E_naive = round(mean(M$naive), 4),
  bias = round(mean(M$naive) - mu, 4))

true_mu E_naive    bias 
 0.2950  0.7769  0.4819

Averaged over 20,000 initial samples the naive mean is 0.7769, a bias of +0.4819, or +163%. That is not a rounding error you can ignore: the naive estimate more than doubles the truth.

The Hansen-Hurwitz estimator

The fix comes from an old idea (Hansen and Hurwitz 1943): weight each observation by how it entered the sample. In adaptive cluster sampling the clean version works through the initial units. A network is included whenever any of its units is drawn initially, and each initial draw points at exactly one network. So assign to each initial unit the mean of its network, and average those network means over the initial sample. Occupied networks contribute their cluster average; empty and edge quadrats contribute their own value (a network of one).

## Hansen-Hurwitz: mean over the initial units of their network mean
hansen_hurwitz <- function(nets, init) mean(nets$net_mean_of_unit[init])
hansen_hurwitz(nets, d1$init)          # the same draw as before

[1] 0.317739

That single draw gives 0.318, close to the truth where the naive mean gave 0.767. The estimator is design-unbiased: because the network means are defined for every unit in the frame and the initial sample is a simple random sample of those units, the Hansen-Hurwitz estimate is just the sample mean of a fixed list of 400 numbers, which is unbiased by construction.

Two overlaid density curves. The naive-mean curve is shifted far to the right of the dashed true-mean line, while the Hansen-Hurwitz curve is centred on the true-mean line. — Figure 3: Sampling distributions over 20,000 initial samples. The naive mean sits well to the right of the truth; the Hansen-Hurwitz estimator is centred on it.

Over the 20,000 samples the Hansen-Hurwitz mean is 0.2947, a bias of -0.00035: unbiased to the precision of the simulation. Its standard deviation is 0.135, against 0.146 for the naive mean, so it is both centred and tighter.

An honest limit on the interval

Unbiasedness is not the whole story. The network means are wildly skewed on a population this sparse, so the normal-based confidence interval built from the standard finite-population variance does not quite reach its nominal cover.

mean(M$cov_hh)     # cover of the nominal 95% Hansen-Hurwitz interval

[1] 0.91005

The interval covers about 91.0% of the time rather than 95%. That is the price of a rare, clumped population: the estimator is honest about the mean but the interval is optimistic, and on a real survey you would lean on a bootstrap or a design with a larger initial sample. The next post keeps the same design but swaps the estimator for the Horvitz-Thompson version, which uses the network inclusion probabilities directly and buys back some of that precision.

References

Hansen MH, Hurwitz WN 1943 Ann Math Stat 14(4):333-362 (10.1214/aoms/1177731356)

Thompson SK 1990 J Am Stat Assoc 85(412):1050-1059 (10.1080/01621459.1990.10474975)

Smith DR, Conroy MJ, Brakhage DH 1995 Biometrics 51(2):777-788 (10.2307/2532964)

Thompson SK, Seber GAF 1996 Adaptive Sampling. Wiley (ISBN 978-0-471-55871-2)

--- title: "Adaptive cluster sampling in R" description: "Build adaptive cluster sampling from scratch in R for rare, clustered species: networks, edge units, why the naive mean is biased, and the Hansen-Hurwitz fix." date: "2026-07-26 14:00" categories: [R, adaptive sampling, survey design, abundance, ecology tutorial] image: thumbnail.png image-alt: "A grid of survey quadrats with a few dense clusters of counts on a mostly empty background, and the adaptive sample growing outward from the occupied cells." --- Rare, clustered species are the hardest thing a fixed design has to cope with. Most quadrats hold nothing, a handful hold almost everything, and a simple random sample of a few dozen units mostly counts zeros. Adaptive cluster sampling turns that structure into an advantage: when a selected unit contains the species, you sample its neighbours too, and keep expanding as long as you keep finding it. The sample follows the organism. The catch is that the mechanism which makes it efficient also breaks the obvious estimator. If you only expand where you find animals, the units you end up observing over-represent the dense patches, so the plain average of everything you sampled runs high. This post builds the design from scratch on a synthetic grid, shows exactly how biased the naive mean is, and derives the Hansen-Hurwitz estimator that puts it right. ```{r} #| label: setup #| include: false options(digits = 6) set.seed(101) ## a rare, clustered count field on a 20x20 grid of quadrats make_population <- function(side = 20L, n_clusters = 5L, per_cluster = 22, spread = 0.9, seed = 101L) { set.seed(seed); N <- side * side; y <- integer(N) cx0 <- runif(n_clusters, 2, side - 1); cy0 <- runif(n_clusters, 2, side - 1) for (k in seq_len(n_clusters)) { np <- rpois(1, per_cluster); if (np == 0) next cx <- pmin(pmax(round(rnorm(np, cx0[k], spread)), 1), side) cy <- pmin(pmax(round(rnorm(np, cy0[k], spread)), 1), side) cell <- (cy - 1) * side + cx tb <- table(cell); idx <- as.integer(names(tb)); y[idx] <- y[idx] + as.integer(tb) } list(side = side, N = N, y = y) } rook <- function(i, side) { x <- ((i - 1) %% side) + 1L; yy <- ((i - 1) %/% side) + 1L nb <- integer(0) if (x > 1) nb <- c(nb, i - 1L) if (x < side) nb <- c(nb, i + 1L) if (yy > 1) nb <- c(nb, i - side) if (yy < side) nb <- c(nb, i + side) nb } build_networks <- function(pop, cval = 1L) { side <- pop$side; N <- pop$N; y <- pop$y; sat <- y >= cval net_id <- integer(N); cur <- 0L; seen <- logical(N) for (i in seq_len(N)) { if (seen[i]) next if (!sat[i]) { cur <- cur + 1L; net_id[i] <- cur; seen[i] <- TRUE; next } cur <- cur + 1L; stack <- i; seen[i] <- TRUE; net_id[i] <- cur while (length(stack)) { u <- stack[length(stack)]; stack <- stack[-length(stack)] for (v in rook(u, side)) if (!seen[v] && sat[v]) { seen[v] <- TRUE; net_id[v] <- cur; stack <- c(stack, v) } } } m <- as.integer(table(net_id)); tot <- as.numeric(tapply(y, net_id, sum)) list(net_id = net_id, m = m, tot = tot, sat = sat, net_mean_of_unit = (tot / m)[net_id], net_is_sat = as.logical(tapply(sat, net_id, function(z) z[1]))) } choose2 <- function(a, b) { out <- exp(lchoose(a, b)); out[b < 0 | b > a] <- 0; out } pop <- make_population() nets <- build_networks(pop, 1L) ubn <- split(seq_len(pop$N), nets$net_id) mu <- mean(pop$y); tau <- sum(pop$y); N <- pop$N occ <- mean(pop$y >= 1) sat_sizes <- sort(nets$m[nets$net_is_sat], decreasing = TRUE) acs_draw <- function(pop, nets, ubn, n1) { N <- pop$N; y <- pop$y; side <- pop$side init <- sample.int(N, n1); hit <- unique(nets$net_id[init]) in_s <- logical(N); in_s[init] <- TRUE for (k in hit) { uk <- ubn[[k]]; in_s[uk] <- TRUE if (nets$net_is_sat[k]) for (u in uk) for (v in rook(u, side)) in_s[v] <- TRUE } a <- vapply(hit, function(k) 1 - choose2(N - nets$m[k], n1) / choose2(N, n1), numeric(1)) list(naive = mean(y[in_s]), hh = mean(nets$net_mean_of_unit[init]), ht = sum(nets$tot[hit] / a) / N, n_final = sum(in_s), init = init, hit = hit) } hh_var <- function(pop, nets, init, n1) { N <- pop$N; w <- nets$net_mean_of_unit[init] (N - n1) / (N * n1 * (n1 - 1)) * sum((w - mean(w))^2) } ``` ## A population built to defeat a fixed design The frame is a grid of `r N` quadrats. The species sits in a few tight clusters on an otherwise empty background: `r sprintf("%.0f%%", 100 * mean(pop$y == 0))` of the quadrats are empty, and the occupancy (the fraction with at least one individual) is only `r sprintf("%.1f%%", 100 * occ)`. The true quadrat mean is `r sprintf("%.3f", mu)` individuals and the true total is `r tau`. On a real survey you would not know these; here they are the yardstick. ```{r} #| label: fig-population #| echo: false #| fig-cap: "The synthetic population. Counts are concentrated in a few clusters on an empty background, the situation adaptive cluster sampling is built for." #| fig-alt: "A 20 by 20 grid of quadrats shaded by count. Most cells are pale (zero), with four small dark clusters of high counts scattered across the grid." #| fig-width: 5.2 #| fig-height: 4.6 library(ggplot2) gd <- data.frame(x = ((seq_len(N) - 1) %% pop$side) + 1, yy = ((seq_len(N) - 1) %/% pop$side) + 1, count = pop$y) ggplot(gd, aes(x, yy, fill = count)) + geom_tile(colour = "#dad9ca", linewidth = 0.25) + scale_fill_gradient(low = "#f5f4ee", high = "#275139", name = "count") + coord_equal() + labs(x = NULL, y = NULL) + theme_minimal(base_size = 11) + theme(panel.grid = element_blank(), axis.text = element_blank(), axis.ticks = element_blank()) ``` Under simple random sampling you would draw a fixed number of quadrats and stop. With `r occ*100`% occupancy most of them are zeros, so the estimate of the mean is dominated by whether you happen to land on a cluster. Adaptive cluster sampling changes the second step: it keeps sampling around any quadrat where the species turns up. ## The design: condition, neighbourhood, networks Three ingredients define the design. The **condition** decides when to expand. Here a quadrat "satisfies the condition" if its count is at least one, that is, the species is present. The **neighbourhood** decides where to expand: we use the four rook neighbours (up, down, left, right). The rule is then recursive. Draw an initial simple random sample; whenever a selected quadrat satisfies the condition, add its neighbours; whenever a newly added quadrat also satisfies the condition, add its neighbours too; stop when no new occupied quadrats appear. The units partition into **networks**. Connected runs of occupied quadrats form one network each: selecting any one of them pulls in all of them. Every quadrat that does not satisfy the condition is a network on its own. The occupied quadrats that sit on the boundary of a network, its non-satisfying neighbours, are the **edge units**: they get observed because they are adjacent to a network, but they do not trigger further expansion. `build_networks()` above finds these by a breadth-first search over the occupied cells. On this population there are `r sum(nets$net_is_sat)` occupied networks, with sizes `r paste(sat_sizes, collapse = ", ")` quadrats. Everything else is a singleton. ```{r} #| label: one-draw #| echo: true ## one adaptive draw from an initial sample of n1 = 30 quadrats set.seed(7) n1 <- 30L d1 <- acs_draw(pop, nets, ubn, n1) c(initial = n1, distinct_networks_hit = length(d1$hit), final_sample = d1$n_final, naive_mean = round(d1$naive, 3)) ``` Thirty initial quadrats expanded to `r d1$n_final` observed quadrats, because a couple of them landed on clusters and dragged the whole cluster (plus its edge) into the sample. The plain average of everything observed in that draw was `r sprintf("%.3f", d1$naive)`, well above the truth of `r sprintf("%.3f", mu)`. ```{r} #| label: fig-expansion #| echo: false #| fig-cap: "One adaptive draw. Open squares are the initial simple random sample; filled squares are quadrats added by expansion. Where an initial unit lands on a cluster, the whole network and its edge come in." #| fig-alt: "The grid with the initial 30 sampled quadrats marked as open squares, scattered at random, and the adaptively added quadrats marked as filled squares, forming compact blocks around two of the clusters." #| fig-width: 5.2 #| fig-height: 4.6 in_s <- logical(N); in_s[d1$init] <- TRUE for (k in d1$hit) { uk <- ubn[[k]]; in_s[uk] <- TRUE if (nets$net_is_sat[k]) for (u in uk) for (v in rook(u, pop$side)) in_s[v] <- TRUE } gd$state <- "not sampled" gd$state[in_s] <- "added by expansion" gd$state[d1$init] <- "initial sample" ggplot(gd, aes(x, yy)) + geom_tile(aes(fill = count), colour = "#dad9ca", linewidth = 0.25) + scale_fill_gradient(low = "#f5f4ee", high = "#cfe0d4", guide = "none") + geom_point(data = subset(gd, state == "added by expansion"), shape = 15, size = 2.2, colour = "#275139") + geom_point(data = subset(gd, state == "initial sample"), shape = 0, size = 2.4, stroke = 0.9, colour = "#b5534e") + coord_equal() + labs(x = NULL, y = NULL) + theme_minimal(base_size = 11) + theme(panel.grid = element_blank(), axis.text = element_blank(), axis.ticks = element_blank()) ``` ## Why the naive mean is biased The bias is not a fluke of one draw. Average the naive mean over many initial samples and it stays high, because expansion is triggered by high values: you always add neighbours around counts, never around empties, so the observed set is enriched for occupied quadrats relative to the frame. Formally the naive mean is not the mean of a probability sample of the frame at all, and it has no reason to be unbiased. ```{r} #| label: mc-naive-hh #| echo: true mc <- function(pop, nets, ubn, n1, B, seed) { set.seed(seed); naive <- hh <- nfin <- numeric(B); cov_hh <- logical(B) mu <- mean(pop$y) for (b in 1:B) { d <- acs_draw(pop, nets, ubn, n1) naive[b] <- d$naive; hh[b] <- d$hh; nfin[b] <- d$n_final se <- sqrt(max(hh_var(pop, nets, d$init, n1), 0)) cov_hh[b] <- abs(d$hh - mu) <= 1.96 * se } list(naive = naive, hh = hh, nfin = nfin, cov_hh = cov_hh) } M <- mc(pop, nets, ubn, n1 = 30L, B = 20000L, seed = 202L) c(true_mu = round(mu, 4), E_naive = round(mean(M$naive), 4), bias = round(mean(M$naive) - mu, 4)) ``` ```{r} #| label: naive-numbers #| include: false naive_bias <- mean(M$naive) - mu naive_pct <- 100 * naive_bias / mu ``` Averaged over 20,000 initial samples the naive mean is `r sprintf("%.4f", mean(M$naive))`, a bias of `r sprintf("%+.4f", naive_bias)`, or `r sprintf("%+.0f%%", naive_pct)`. That is not a rounding error you can ignore: the naive estimate more than doubles the truth. ## The Hansen-Hurwitz estimator The fix comes from an old idea (Hansen and Hurwitz 1943): weight each observation by how it entered the sample. In adaptive cluster sampling the clean version works through the initial units. A network is included whenever any of its units is drawn initially, and each initial draw points at exactly one network. So assign to each initial unit the **mean of its network**, and average those network means over the initial sample. Occupied networks contribute their cluster average; empty and edge quadrats contribute their own value (a network of one). ```{r} #| label: hh-definition #| echo: true ## Hansen-Hurwitz: mean over the initial units of their network mean hansen_hurwitz <- function(nets, init) mean(nets$net_mean_of_unit[init]) hansen_hurwitz(nets, d1$init) # the same draw as before ``` That single draw gives `r sprintf("%.3f", d1$hh)`, close to the truth where the naive mean gave `r sprintf("%.3f", d1$naive)`. The estimator is design-unbiased: because the network means are defined for every unit in the frame and the initial sample is a simple random sample of those units, the Hansen-Hurwitz estimate is just the sample mean of a fixed list of `r N` numbers, which is unbiased by construction. ```{r} #| label: hh-checks #| include: false hh_bias <- mean(M$hh) - mu ``` ```{r} #| label: fig-naive-vs-hh #| echo: false #| fig-cap: "Sampling distributions over 20,000 initial samples. The naive mean sits well to the right of the truth; the Hansen-Hurwitz estimator is centred on it." #| fig-alt: "Two overlaid density curves. The naive-mean curve is shifted far to the right of the dashed true-mean line, while the Hansen-Hurwitz curve is centred on the true-mean line." #| fig-width: 6.0 #| fig-height: 3.7 dd <- rbind(data.frame(est = M$naive, which = "naive mean"), data.frame(est = M$hh, which = "Hansen-Hurwitz")) ggplot(dd, aes(est, colour = which)) + geom_density(linewidth = 0.9) + geom_vline(xintercept = mu, linetype = 2, colour = "#16241d") + scale_colour_manual(values = c("naive mean" = "#b5534e", "Hansen-Hurwitz" = "#275139"), name = NULL) + labs(x = "estimated mean count per quadrat", y = "density") + coord_cartesian(xlim = c(0, 1.3)) + theme_minimal(base_size = 11) + theme(legend.position = "top", panel.grid.minor = element_blank()) ``` Over the 20,000 samples the Hansen-Hurwitz mean is `r sprintf("%.4f", mean(M$hh))`, a bias of `r sprintf("%+.5f", hh_bias)`: unbiased to the precision of the simulation. Its standard deviation is `r sprintf("%.3f", sd(M$hh))`, against `r sprintf("%.3f", sd(M$naive))` for the naive mean, so it is both centred and tighter. ## An honest limit on the interval Unbiasedness is not the whole story. The network means are wildly skewed on a population this sparse, so the normal-based confidence interval built from the standard finite-population variance does not quite reach its nominal cover. ```{r} #| label: hh-coverage #| echo: true mean(M$cov_hh) # cover of the nominal 95% Hansen-Hurwitz interval ``` The interval covers about `r sprintf("%.1f%%", 100 * mean(M$cov_hh))` of the time rather than 95%. That is the price of a rare, clumped population: the estimator is honest about the mean but the interval is optimistic, and on a real survey you would lean on a bootstrap or a design with a larger initial sample. The next post keeps the same design but swaps the estimator for the Horvitz-Thompson version, which uses the network inclusion probabilities directly and buys back some of that precision. ## Related tutorials - [Horvitz-Thompson for adaptive samples](../horvitz-thompson-for-adaptive-samples/) - [When does adaptive sampling pay?](../when-does-adaptive-sampling-pay/) - [Checking an adaptive sampling design](../checking-an-adaptive-sampling-design/) - [Stratified random sampling in ecology](../stratified-random-sampling/) ## References Hansen MH, Hurwitz WN 1943 Ann Math Stat 14(4):333-362 (10.1214/aoms/1177731356) Thompson SK 1990 J Am Stat Assoc 85(412):1050-1059 (10.1080/01621459.1990.10474975) Smith DR, Conroy MJ, Brakhage DH 1995 Biometrics 51(2):777-788 (10.2307/2532964) Thompson SK, Seber GAF 1996 Adaptive Sampling. Wiley (ISBN 978-0-471-55871-2)