Horvitz-Thompson for adaptive samples

Estimate abundance from an adaptive cluster sample with the Horvitz-Thompson estimator: network inclusion probabilities in R, and why it beats Hansen-Hurwitz.

Corrected 14 October 2026: the post said neither estimator is a function of the minimal sufficient statistic, so that the Rao-Blackwell step lowers the variance of both; on this design, where the sample expands from every quadrat holding the species, that holds for Hansen-Hurwitz only, and the post now says the step leaves Horvitz-Thompson unchanged. It also no longer says that field guidance leans on a bootstrap interval, a claim the post did not support.

The previous post built adaptive cluster sampling and fixed the biased naive mean with the Hansen-Hurwitz estimator. Hansen-Hurwitz is unbiased and easy to compute, but it is not the tightest estimator the design allows. This post derives the alternative: the Horvitz-Thompson estimator, which works from the probability that each network makes it into the sample. It is a little more work to set up and it wins on precision.

The idea is the same one that runs through all unequal-probability sampling (Horvitz and Thompson 1952): if you can write down the probability that each thing you observed had of being observed, you can weight each observation by the inverse of that probability and recover an unbiased total. In adaptive cluster sampling the “things” are networks, and their inclusion probabilities have a clean closed form.

options(digits = 6)

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); N <- pop$N; n1 <- 30L
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) {
  w <- nets$net_mean_of_unit[init]
  (N - n1) / (N * n1 * (n1 - 1)) * sum((w - mean(w))^2)
}
ht_var <- function(pop, nets, hit, n1) {
  K <- length(hit); if (K < 2) return(NA_real_)
  m <- nets$m[hit]; tot <- nets$tot[hit]; Ctot <- choose2(N, n1)
  a <- 1 - vapply(m, function(mk) choose2(N - mk, n1), numeric(1)) / Ctot
  v <- 0
  for (j in 1:K) for (l in 1:K) {
    ajl <- if (j == l) a[j] else
      1 - (choose2(N - m[j], n1) + choose2(N - m[l], n1) -
             choose2(N - m[j] - m[l], n1)) / Ctot
    v <- v + (ajl - a[j] * a[l]) / ajl * (tot[j] / a[j]) * (tot[l] / a[l])
  }
  v / N^2
}

The probability a network gets into the sample

A network of m units is included in the sample if and only if the initial simple random sample lands on at least one of its units. The complement is easier to count: the probability the initial sample of 30 units avoids the whole network. So the inclusion probability of a network of size \(m\) is

\[ \alpha = 1 - \frac{\binom{N-m}{n_1}}{\binom{N}{n_1}}, \]

with \(N\) the number of quadrats and \(n_1\) the initial sample size. Bigger networks are far more likely to be hit; a single quadrat is included only with the plain initial probability \(n_1/N\).

alpha <- function(m, N, n1) 1 - choose2(N - m, n1) / choose2(N, n1)
m_big <- max(nets$m)                        # the largest occupied network
c(largest_network = m_big,
  alpha_largest = round(alpha(m_big, N, n1), 4),
  alpha_single  = round(alpha(1, N, n1), 4),   # = n1 / N
  n1_over_N     = round(n1 / N, 4))

largest_network   alpha_largest    alpha_single       n1_over_N 
        19.0000          0.7806          0.0750          0.0750 

The largest occupied network here spans 19 quadrats and has inclusion probability 0.7806: with 30 initial draws it is caught almost four times out of five. A lone quadrat has probability only 0.0750, exactly \(n_1/N\). That spread in inclusion probabilities is what the estimator has to correct for.

Figure 1: Inclusion probability against network size for the initial sample of 30 quadrats. A single quadrat enters with probability n1/N; large clusters are nearly certain to be caught.

The Horvitz-Thompson estimator

With inclusion probabilities in hand, the estimator writes itself. Take the distinct networks the initial sample intersected, divide each network’s total by its inclusion probability, add them up, and divide by the number of quadrats for the mean.

horvitz_thompson <- function(nets, hit, N, n1) {
  a <- alpha(nets$m[hit], N, n1)
  sum(nets$tot[hit] / a) / N
}
set.seed(7); d1 <- acs_draw(pop, nets, ubn, n1)
horvitz_thompson(nets, d1$hit, N, n1)

[1] 0.244039

The word distinct is doing real work. If two of your initial units land in the same cluster, that network was still included once, so it enters the sum once, weighted by its single inclusion probability. This is where Horvitz-Thompson parts company with Hansen-Hurwitz, which effectively counts a network once per initial hit.

Why it beats Hansen-Hurwitz

On a clustered population the initial sample often hits the same network more than once: there are only a few networks and 30 initial draws.

set.seed(202); B <- 20000L; multi <- logical(B)
for (b in 1:B) {
  init <- sample.int(N, n1); tab <- table(nets$net_id[init])
  multi[b] <- any(tab > 1 & nets$net_is_sat[as.integer(names(tab))])
}
mean(multi)          # fraction of draws with a multiply-hit occupied network

[1] 0.71445

An occupied network is hit more than once in 71% of draws. Hansen-Hurwitz counts each of those hits, so the same cluster enters the average twice, and sometimes more often than that, and the repetition adds variance for no information. Horvitz-Thompson counts the cluster once. In estimation-theory terms, Hansen-Hurwitz is not a function of the minimal sufficient statistic. That statistic is the unordered set of distinct quadrats observed together with their counts, and Hansen-Hurwitz additionally needs to know which quadrats were in the initial sample, something the final data do not record. Averaging it over all the initial samples consistent with the observed data, the Rao-Blackwell step, therefore lowers its variance (Salehi 1999). Horvitz-Thompson needs no such record on this design: the sample expands from every quadrat where the species is present, so every quadrat outside an occupied network holds zero, the only networks that add to its sum are the occupied networks in the final sample, and the same step leaves it unchanged. Under a stricter rule for expanding, a quadrat that does not trigger expansion can still hold a nonzero count, and then Horvitz-Thompson too depends on whether that quadrat was in the initial sample.

mc <- function(pop, nets, ubn, n1, B, seed) {
  set.seed(seed); hh <- ht <- numeric(B); cov_ht <- logical(B); neg <- 0L
  mu <- mean(pop$y)
  for (b in 1:B) {
    d <- acs_draw(pop, nets, ubn, n1)
    hh[b] <- d$hh; ht[b] <- d$ht
    vt <- ht_var(pop, nets, d$hit, n1)
    if (!is.na(vt) && vt < 0) neg <- neg + 1L
    cov_ht[b] <- if (is.na(vt) || vt <= 0) NA else abs(d$ht - mu) <= 1.96 * sqrt(vt)
  }
  list(hh = hh, ht = ht, cov_ht = cov_ht, neg = neg)
}
M <- mc(pop, nets, ubn, n1 = 30L, B = 20000L, seed = 202L)
c(bias_ht = round(mean(M$ht) - mu, 5),
  sd_hh = round(sd(M$hh), 4), sd_ht = round(sd(M$ht), 4),
  var_ratio_ht_over_hh = round(var(M$ht) / var(M$hh), 3))

             bias_ht                sd_hh                sd_ht 
            -0.00003              0.13490              0.10370 
var_ratio_ht_over_hh 
             0.59100 

ratio <- var(M$ht) / var(M$hh)
se_ratio <- sd(M$ht) / sd(M$hh)
ht_cover <- mean(M$cov_ht, na.rm = TRUE)

Both estimators are unbiased: the Horvitz-Thompson bias over 20,000 samples is -0.00003. But the Horvitz-Thompson variance is 0.591 times the Hansen-Hurwitz variance, a standard error of 0.1037 against 0.1349, or about 23% narrower. On this population the improvement is worth having, and Salehi (2003) recommends Horvitz-Thompson as the default despite the heavier bookkeeping.

Figure 2: Sampling distributions of the two unbiased estimators over 20,000 initial samples. Both centre on the truth; the Horvitz-Thompson estimator is the tighter of the two.

The interval is still the weak point

More precise does not mean better calibrated. The unbiased Horvitz-Thompson variance estimator uses pairwise network inclusion probabilities, and on a skewed population the normal interval built from it undercovers, a touch worse than the Hansen-Hurwitz one from the previous post.

c(ht_interval_cover = round(ht_cover, 3),
  negative_variance_estimates = M$neg)

          ht_interval_cover negative_variance_estimates 
                      0.883                       0.000 

The interval reaches only about 88.3% cover here, and that average is taken only over the draws with a positive variance estimate; in the remaining 1.3% of draws the initial sample missed every occupied network, so the estimate and its estimated variance are both zero, the interval collapses to a point, and those draws are dropped rather than scored as misses. No variance estimate came out negative on this population, but the Horvitz-Thompson variance estimator can go negative in general, so a survey that relies on it has to check for that. The point estimate is excellent; the interval needs care. The next post steps back from the estimator to the design itself and asks the question that decides whether any of this is worth doing: when does adaptive sampling actually beat a simple random sample?

References

Horvitz DG, Thompson DJ 1952 J Am Stat Assoc 47(260):663-685 (10.1080/01621459.1952.10483446)

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

Salehi MM 1999 Environ Ecol Stat 6(2):183-195 (10.1023/A:1009670205509)

Salehi MM 2003 Environ Ecol Stat 10(1):115-127 (10.1023/A:1021989509323)

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

Related tutorials

• Adaptive cluster sampling in R
• When does adaptive sampling pay?
• Checking an adaptive sampling design
• Unequal-probability spatial sampling

Reproducibility

• Code last run: 2026-10-15
• R version 4.5.3 (2026-03-11), aarch64-apple-darwin20
• Packages: ggplot2 4.0.3