Two-phase sampling when strata are unknown

survey design

sampling

ecology tutorial

Double sampling in R: a large cheap first phase builds the strata, a small expensive second phase measures the response, and the budget split follows a formula.

Author

Tidy Ecology

Published

2026-07-26

Stratified sampling assumes a map. You need to know which stratum every plot belongs to and how large each stratum is, before a single measurement is taken. Plenty of surveys have no such map: the habitat classification is exactly what you would have to go out and produce, and by then you have spent the season.

Two-phase sampling, also called double sampling, buys the map with part of the budget. Phase one visits a large sample and records something cheap: a habitat call from the track, a canopy class, a rapid cover estimate. Phase two returns to a subsample of those units and does the expensive thing: harvesting, weighing, sorting, sequencing. The cheap phase supplies the stratum weights, the expensive phase supplies the stratum means, and the two are combined into one estimate.

The interesting part is that the variance splits neatly along the same seam, which turns the budget split into arithmetic instead of a hunch.

The setup

Same landscape, but now nobody has the habitat map:

set.seed(371)

Nh   <- c(meadow = 600, wetland = 300, woodland = 300)
mu_h <- c(meadow = 12,  wetland = 26,  woodland = 5)
sd_h <- c(meadow = 4,   wetland = 9,   woodland = 2)

N       <- sum(Nh)
stratum <- rep(names(Nh), Nh)
y       <- round(pmax(0, rnorm(N, rep(mu_h, Nh), rep(sd_h, Nh))), 2)

W   <- Nh / N
mu  <- mean(y)
S2  <- var(y)
S2h <- tapply(y, stratum, var)[names(Nh)]

within  <- sum(W * S2h)          # what stays inside the strata
between <- S2 - within           # what the strata remove
c(total = S2, between = between, within = within)

   total  between   within 
85.92611 58.08970 27.83641

68 percent of the variance sits between the strata and the rest sits inside them. Hold on to those two numbers, because they turn out to be the price tags of the two phases.

The estimator

Draw a simple random sample of n1 units and record the cheap classification. That gives estimated weights $w_h = n_{1h}/n_1$. Then subsample each phase-one stratum and measure the response, giving $\bar y_h$. The estimator is the same weighted sum as before, with estimated weights in place of known ones:

\[\bar y_{\text{2ph}} = \sum_h w_h \bar y_h .\]

It is unbiased. Conditional on the phase-one sample, the expectation is the phase-one sample mean of the response, and the phase-one sample mean is unbiased for the population mean, so the two steps compose.

The variance splits along the phases

Condition on phase one and add the two pieces. The variation of the phase-one mean contributes $(1 - n_1/N)S^2/n_1$; the subsampling contributes the usual within-stratum term, evaluated at the phase-two sizes. Written with a common subsampling fraction it collapses to

\[V(\bar y_{\text{2ph}}) \;\approx\; \Big(1 - \tfrac{n_1}{N}\Big)\frac{S^2}{n_1} \;+\; \Big(\sum_h W_h S_h^2\Big)\left(\frac{1}{n_2} - \frac{1}{n_1}\right),\]

and, ignoring the finite population correction, into something worth putting on a wall:

\[V \;\approx\; \frac{\text{between}}{n_1} \;+\; \frac{\text{within}}{n_2}.\]

The cheap phase pays for the between-stratum variance and the expensive phase pays for the within-stratum variance. That is the whole design in one line.

V_2ph <- function(n1, n2) (1 - n1 / N) * S2 / n1 + within * (1 / n2 - 1 / n1)
c(`n1=400, n2=100` = V_2ph(400, 100), `n1=300, n2=60` = V_2ph(300, 60))

n1=400, n2=100  n1=300, n2=60 
     0.3519833      0.5859674

Does the formula hold?

Simulate the whole two-phase procedure, weights and all, and compare.

one_survey <- function(n1, n2) {
  s1 <- sample.int(N, n1)
  st <- stratum[s1]
  w  <- table(factor(st, levels = names(Nh))) / n1
  nu <- n2 / n1
  est <- 0
  for (h in names(Nh)) {
    ix <- s1[st == h]
    est <- est + w[[h]] * mean(y[sample(ix, max(1, round(nu * length(ix))))])
  }
  est
}

B <- 3000
res <- t(sapply(list(c(400, 100), c(300, 60)), function(p) {
  set.seed(3731)
  mc <- replicate(B, one_survey(p[1], p[2]))
  c(n1 = p[1], n2 = p[2], mc_mean = mean(mc), truth = mu,
    mc_variance = var(mc), formula = V_2ph(p[1], p[2]))
}))
round(res, 4)

      n1  n2 mc_mean   truth mc_variance formula
[1,] 400 100 13.8164 13.7966      0.3395   0.352
[2,] 300  60 13.8339 13.7966      0.5638   0.586

The simulated variance is 0.340 against a predicted 0.352, so the formula runs about 4 percent high. That direction is expected: the derivation replaces the random phase-one stratum counts with their expectations, and the approximation is a large-sample one. The estimator itself is centred, with a simulated bias of 0.020 on a mean of 13.80.

Splitting a budget

Let phase one cost $c_1$ per unit and phase two cost $c_2$, and fix the total. Minimising $A/n_1 + B/n_2$ under $c_1 n_1 + c_2 n_2 = C$ is a one-line Lagrange problem with a memorable answer,

\[\frac{n_1}{n_2} = \sqrt{\frac{A\,c_2}{B\,c_1}},\]

where $A$ is the between-stratum component and $B$ the within-stratum one. Spend more on the phase that is buying the larger piece of variance, discounted by what that phase costs.

c1 <- 1; c2 <- 10; budget <- 1000

ratio <- sqrt(between * c2 / (within * c1))
n2_opt <- budget / (ratio * c1 + c2)
n1_opt <- ratio * n2_opt
c(n1 = n1_opt, n2 = n2_opt, V = V_2ph(n1_opt, n2_opt), SE = sqrt(V_2ph(n1_opt, n2_opt)))

         n1          n2           V          SE 
313.5723397  68.6427660   0.5191720   0.7205359

U-shaped curve of design variance against the share of budget spent on phase one, with a minimum around one third and a marked optimum point. — Figure 1: Design variance against how the budget is divided, holding total spend fixed. Too little phase one and the weights are noisy; too much and there is nothing left to measure the response with.

What the design is worth

Three ways to spend the same money. Spend it all on the expensive measurement and take a simple random sample. Split it between the phases. Or, the unattainable benchmark, have the map already and run a proper stratified survey.

n_srs <- budget / c2
V_srs <- (1 - n_srs / N) * S2 / n_srs
V_map <- (1 - n_srs / N) / n_srs * within        # stratified, weights known in advance

comparison <- rbind(
  `simple random`  = c(expensive_units = n_srs, V = V_srs,                SE = sqrt(V_srs)),
  `two-phase`      = c(expensive_units = n2_opt, V = V_2ph(n1_opt, n2_opt),
                       SE = sqrt(V_2ph(n1_opt, n2_opt))),
  `map already in hand` = c(expensive_units = n_srs, V = V_map,           SE = sqrt(V_map)))
round(comparison, 4)

                    expensive_units      V     SE
simple random              100.0000 0.7877 0.8875
two-phase                   68.6428 0.5192 0.7205
map already in hand        100.0000 0.2552 0.5051

Three bars of standard error: simple random highest, two-phase in the middle, map-in-hand stratified lowest. — Figure 2: Standard error on an identical budget. Two-phase sampling recovers part of the distance between a simple random sample and a survey that already had the map, and it recovers it without a map.

Two-phase sampling closes 50 percent of the variance gap between having no map and having one, standard error 0.721 against 0.887 and 0.505. Buying the map with survey effort is worth doing, and it is not the same as being handed one.

When the cheap phase is not cheap enough

All of that depends on $c_1$. Push the first phase towards the price of the second and the design stops paying for itself, because every phase-one unit is a phase-two unit you did not buy.

V_at <- function(cc) {
  r <- sqrt(between * c2 / (within * cc))
  n2 <- budget / (r * cc + c2)
  V_2ph(r * n2, n2)
}
break_even <- uniroot(function(cc) V_at(cc) - V_srs, c(0.2, 20))$root

sweep <- data.frame(c1 = c(0.2, 0.5, 1, 2, 4))
sweep$V     <- vapply(sweep$c1, V_at, 0)
sweep$ratio <- sweep$V / V_srs
round(sweep, 4)

   c1      V  ratio
1 0.2 0.3321 0.4216
2 0.5 0.4156 0.5277
3 1.0 0.5192 0.6591
4 2.0 0.6826 0.8666
5 4.0 0.9478 1.2033

c(break_even_c1 = break_even, as_fraction_of_c2 = break_even / c2)

    break_even_c1 as_fraction_of_c2 
        2.7455658         0.2745566

The break-even point is a phase-one unit costing 27 percent of a phase-two unit. Below that, two phases beat one; above it, you should have spent everything on the response. In this landscape the classification would have to be genuinely quick, a call from the vehicle rather than a walked transect.

Rising curve of relative variance against phase-one cost, crossing the value one just before a phase-one cost of three. — Figure 3: Two-phase variance relative to a simple random sample, as the cost of a phase-one unit rises. The design pays for itself only to the left of the crossing point.

Two refinements worth knowing

The subsampling fraction does not have to be the same in every stratum. Minimising the same variance with $\nu_h$ free gives the Neyman analogue, $\nu_h \propto S_h$, and the within term becomes $(\sum_h W_h S_h)^2/n_2$:

Sh <- sqrt(S2h)
V_2ph_opt <- function(n1, n2) (1 - n1 / N) * S2 / n1 + (sum(W * Sh))^2 / n2 - within / n1
r2  <- sqrt((S2 - (sum(W * Sh))^2) * c2 / ((sum(W * Sh))^2 * c1))
n2b <- budget / (r2 * c1 + c2); n1b <- r2 * n2b
c(n1 = n1b, n2 = n2b, V_uniform = V_2ph(n1_opt, n2_opt), V_optimal_nu = V_2ph_opt(n1b, n2b))

          n1           n2    V_uniform V_optimal_nu 
 352.8117022   64.7188298    0.5191720    0.4273174

Subsampling in proportion to the stratum standard deviations takes another 18 percent off the variance, for free, using information the first phase has already collected. It is the same Neyman argument one level down.

The second refinement is a warning. The phase-one classification is allowed to be wrong. Because the estimator weights whatever classes the field call produced, misclassification does not bias anything; it only makes the classes separate the response less well, which shows up as a smaller between component.

set.seed(3733)
proxy <- do.call(rbind, lapply(c(0, 0.10, 0.20, 0.35), function(e) {
  cls <- stratum
  flip <- runif(N) < e
  cls[flip] <- sample(names(Nh), sum(flip), replace = TRUE)
  Wc  <- as.numeric(table(factor(cls, levels = names(Nh)))) / N
  wi  <- sum(Wc * tapply(y, cls, var)[names(Nh)])
  r   <- sqrt((S2 - wi) * c2 / (wi * c1)); n2e <- budget / (r * c1 + c2); n1e <- r * n2e
  Ve  <- (1 - n1e / N) * S2 / n1e + wi * (1 / n2e - 1 / n1e)
  data.frame(error = e, within_share = wi / S2, V = Ve, ratio_to_srs = Ve / V_srs)
}))
round(proxy, 3)

  error within_share     V ratio_to_srs
1  0.00        0.324 0.519        0.659
2  0.10        0.483 0.659        0.837
3  0.20        0.551 0.711        0.903
4  0.35        0.745 0.827        1.051

A tenth of the units misclassified moves the relative variance from 0.66 to 0.84, and at 35 percent misclassification the design has slipped past parity to 1.05. Validity is safe, precision is not. A cheap classification that is cheap because it is careless will quietly spend your budget on nothing.

The honest limit

Everything here rests on phase one being a probability sample of the same frame as phase two, with phase two drawn inside it. Convenience units bolted on at either stage break the chain, and no amount of careful weighting repairs them.

The variance formula is an approximation valid when the phase-one stratum counts are comfortably large. With a small first phase, or a class so rare that it collects a handful of units, the approximation runs optimistic in the wrong direction and the realised precision is worse than advertised. Two-phase designs are also more work to run than they look on paper: revisiting a specific subsample means keeping track of which units were which, in the field, months apart.

Where to go next

Three posts of design and no diagnostics yet. The last post in this group takes a finished stratified survey and asks the awkward questions: what happens when you stratify after the fact instead of before, what an out-of-date area figure does to the estimate, and how many strata are too many.

References

Neyman 1938 J Am Stat Assoc 33(201):101-116 (10.1080/01621459.1938.10503378)

Rao 1973 Biometrika 60(1):125-133 (10.1093/biomet/60.1.125)

Cochran 1977 Sampling Techniques, 3rd edn, Wiley; ISBN 978-0-471-16240-7

Sarndal, Swensson and Wretman 1992 Model Assisted Survey Sampling, Springer; ISBN 978-0-387-40620-6

Gregoire and Valentine 2008 Sampling Strategies for Natural Resources and the Environment, Chapman and Hall; ISBN 978-1-58488-370-8

--- title: "Two-phase sampling when strata are unknown" description: "Double sampling in R: a large cheap first phase builds the strata, a small expensive second phase measures the response, and the budget split follows a formula." date: "2026-07-26 12:00" categories: [R, survey design, sampling, ecology tutorial] image: thumbnail.png image-alt: "U-shaped curve of design variance against the share of budget spent on the cheap first phase, with the optimum marked" --- Stratified sampling assumes a map. You need to know which stratum every plot belongs to and how large each stratum is, before a single measurement is taken. Plenty of surveys have no such map: the habitat classification is exactly what you would have to go out and produce, and by then you have spent the season. Two-phase sampling, also called double sampling, buys the map with part of the budget. Phase one visits a large sample and records something cheap: a habitat call from the track, a canopy class, a rapid cover estimate. Phase two returns to a subsample of those units and does the expensive thing: harvesting, weighing, sorting, sequencing. The cheap phase supplies the stratum weights, the expensive phase supplies the stratum means, and the two are combined into one estimate. The interesting part is that the variance splits neatly along the same seam, which turns the budget split into arithmetic instead of a hunch. ## The setup Same landscape, but now nobody has the habitat map: ```{r frame} set.seed(371) Nh <- c(meadow = 600, wetland = 300, woodland = 300) mu_h <- c(meadow = 12, wetland = 26, woodland = 5) sd_h <- c(meadow = 4, wetland = 9, woodland = 2) N <- sum(Nh) stratum <- rep(names(Nh), Nh) y <- round(pmax(0, rnorm(N, rep(mu_h, Nh), rep(sd_h, Nh))), 2) W <- Nh / N mu <- mean(y) S2 <- var(y) S2h <- tapply(y, stratum, var)[names(Nh)] within <- sum(W * S2h) # what stays inside the strata between <- S2 - within # what the strata remove c(total = S2, between = between, within = within) ``` ```{r scalars-0} #| echo: false share_between <- between / S2 ``` `r sprintf("%.0f", 100 * share_between)` percent of the variance sits between the strata and the rest sits inside them. Hold on to those two numbers, because they turn out to be the price tags of the two phases. ## The estimator Draw a simple random sample of `n1` units and record the cheap classification. That gives estimated weights $w_h = n_{1h}/n_1$. Then subsample each phase-one stratum and measure the response, giving $\bar y_h$. The estimator is the same weighted sum as before, with estimated weights in place of known ones: $$\bar y_{\text{2ph}} = \sum_h w_h \bar y_h .$$ It is unbiased. Conditional on the phase-one sample, the expectation is the phase-one sample mean of the response, and the phase-one sample mean is unbiased for the population mean, so the two steps compose. ## The variance splits along the phases Condition on phase one and add the two pieces. The variation of the phase-one mean contributes $(1 - n_1/N)S^2/n_1$; the subsampling contributes the usual within-stratum term, evaluated at the phase-two sizes. Written with a common subsampling fraction it collapses to $$V(\bar y_{\text{2ph}}) \;\approx\; \Big(1 - \tfrac{n_1}{N}\Big)\frac{S^2}{n_1} \;+\; \Big(\sum_h W_h S_h^2\Big)\left(\frac{1}{n_2} - \frac{1}{n_1}\right),$$ and, ignoring the finite population correction, into something worth putting on a wall: $$V \;\approx\; \frac{\text{between}}{n_1} \;+\; \frac{\text{within}}{n_2}.$$ **The cheap phase pays for the between-stratum variance and the expensive phase pays for the within-stratum variance.** That is the whole design in one line. ```{r variance-fn} V_2ph <- function(n1, n2) (1 - n1 / N) * S2 / n1 + within * (1 / n2 - 1 / n1) c(`n1=400, n2=100` = V_2ph(400, 100), `n1=300, n2=60` = V_2ph(300, 60)) ``` ## Does the formula hold? Simulate the whole two-phase procedure, weights and all, and compare. ```{r monte-carlo} one_survey <- function(n1, n2) { s1 <- sample.int(N, n1) st <- stratum[s1] w <- table(factor(st, levels = names(Nh))) / n1 nu <- n2 / n1 est <- 0 for (h in names(Nh)) { ix <- s1[st == h] est <- est + w[[h]] * mean(y[sample(ix, max(1, round(nu * length(ix))))]) } est } B <- 3000 res <- t(sapply(list(c(400, 100), c(300, 60)), function(p) { set.seed(3731) mc <- replicate(B, one_survey(p[1], p[2])) c(n1 = p[1], n2 = p[2], mc_mean = mean(mc), truth = mu, mc_variance = var(mc), formula = V_2ph(p[1], p[2])) })) round(res, 4) ``` ```{r scalars-1} #| echo: false mc_v1 <- res[1, "mc_variance"]; fo_v1 <- res[1, "formula"] mc_b1 <- res[1, "mc_mean"] - mu rel1 <- 100 * (fo_v1 / mc_v1 - 1) ``` The simulated variance is `r sprintf("%.3f", mc_v1)` against a predicted `r sprintf("%.3f", fo_v1)`, so the formula runs about `r sprintf("%.0f", rel1)` percent high. That direction is expected: the derivation replaces the random phase-one stratum counts with their expectations, and the approximation is a large-sample one. The estimator itself is centred, with a simulated bias of `r sprintf("%.3f", mc_b1)` on a mean of `r sprintf("%.2f", mu)`. ## Splitting a budget Let phase one cost $c_1$ per unit and phase two cost $c_2$, and fix the total. Minimising $A/n_1 + B/n_2$ under $c_1 n_1 + c_2 n_2 = C$ is a one-line Lagrange problem with a memorable answer, $$\frac{n_1}{n_2} = \sqrt{\frac{A\,c_2}{B\,c_1}},$$ where $A$ is the between-stratum component and $B$ the within-stratum one. Spend more on the phase that is buying the larger piece of variance, discounted by what that phase costs. ```{r budget-split} c1 <- 1; c2 <- 10; budget <- 1000 ratio <- sqrt(between * c2 / (within * c1)) n2_opt <- budget / (ratio * c1 + c2) n1_opt <- ratio * n2_opt c(n1 = n1_opt, n2 = n2_opt, V = V_2ph(n1_opt, n2_opt), SE = sqrt(V_2ph(n1_opt, n2_opt))) ``` ```{r fig-split} #| echo: false #| fig-cap: "Design variance against how the budget is divided, holding total spend fixed. Too little phase one and the weights are noisy; too much and there is nothing left to measure the response with." #| fig-alt: "U-shaped curve of design variance against the share of budget spent on phase one, with a minimum around one third and a marked optimum point." library(ggplot2) theme_te <- function() { theme_minimal(base_size = 11) + theme( panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#dad9ca", linewidth = 0.3), plot.title = element_text(colour = "#16241d", face = "bold", size = 12), plot.subtitle = element_text(colour = "#5d6b61", size = 9.5), axis.title = element_text(colour = "#46604a", size = 9.5), axis.text = element_text(colour = "#2c3a31", size = 9), legend.title = element_text(colour = "#46604a", size = 9.5), legend.text = element_text(colour = "#2c3a31", size = 9), strip.text = element_text(colour = "#16241d", face = "bold", size = 10) ) } sh <- seq(0.05, 0.85, by = 0.01) cv <- data.frame(share = sh, V = vapply(sh, function(s) V_2ph(budget * s / c1, budget * (1 - s) / c2), 0)) ggplot(cv, aes(share, V)) + geom_line(colour = "#275139", linewidth = 0.9) + annotate("point", x = n1_opt * c1 / budget, y = V_2ph(n1_opt, n2_opt), colour = "#b5534e", size = 3) + annotate("text", x = n1_opt * c1 / budget, y = V_2ph(n1_opt, n2_opt), label = " optimum", hjust = 0, vjust = 2, colour = "#b5534e", size = 3.2) + labs(title = "The budget split has an interior optimum", subtitle = "Phase one costs a tenth of phase two; total spend held fixed", x = "share of the budget spent on phase one", y = "design variance of the mean") + theme_te() ``` ## What the design is worth Three ways to spend the same money. Spend it all on the expensive measurement and take a simple random sample. Split it between the phases. Or, the unattainable benchmark, have the map already and run a proper stratified survey. ```{r three-way} n_srs <- budget / c2 V_srs <- (1 - n_srs / N) * S2 / n_srs V_map <- (1 - n_srs / N) / n_srs * within # stratified, weights known in advance comparison <- rbind( `simple random` = c(expensive_units = n_srs, V = V_srs, SE = sqrt(V_srs)), `two-phase` = c(expensive_units = n2_opt, V = V_2ph(n1_opt, n2_opt), SE = sqrt(V_2ph(n1_opt, n2_opt))), `map already in hand` = c(expensive_units = n_srs, V = V_map, SE = sqrt(V_map))) round(comparison, 4) ``` ```{r scalars-2} #| echo: false gap_closed <- (V_srs - V_2ph(n1_opt, n2_opt)) / (V_srs - V_map) se_srs <- sqrt(V_srs); se_2ph <- sqrt(V_2ph(n1_opt, n2_opt)); se_map <- sqrt(V_map) ``` ```{r fig-three} #| echo: false #| fig-cap: "Standard error on an identical budget. Two-phase sampling recovers part of the distance between a simple random sample and a survey that already had the map, and it recovers it without a map." #| fig-alt: "Three bars of standard error: simple random highest, two-phase in the middle, map-in-hand stratified lowest." cd <- data.frame(design = factor(rownames(comparison), levels = rownames(comparison)), SE = comparison[, "SE"]) ggplot(cd, aes(design, SE)) + geom_col(fill = c("#b5534e", "#93a87f", "#275139"), width = 0.55) + geom_text(aes(label = sprintf("%.3f", SE)), vjust = -0.6, colour = "#5d6b61", size = 3.2) + scale_y_continuous(expand = expansion(mult = c(0, 0.14))) + labs(title = "What the cheap phase buys you", subtitle = "Same field budget in all three; only the middle bar is achievable without a map", x = NULL, y = "standard error of the mean") + theme_te() ``` Two-phase sampling closes `r sprintf("%.0f", 100 * gap_closed)` percent of the variance gap between having no map and having one, standard error `r sprintf("%.3f", se_2ph)` against `r sprintf("%.3f", se_srs)` and `r sprintf("%.3f", se_map)`. Buying the map with survey effort is worth doing, and it is not the same as being handed one. ## When the cheap phase is not cheap enough All of that depends on $c_1$. Push the first phase towards the price of the second and the design stops paying for itself, because every phase-one unit is a phase-two unit you did not buy. ```{r cost-sweep} V_at <- function(cc) { r <- sqrt(between * c2 / (within * cc)) n2 <- budget / (r * cc + c2) V_2ph(r * n2, n2) } break_even <- uniroot(function(cc) V_at(cc) - V_srs, c(0.2, 20))$root sweep <- data.frame(c1 = c(0.2, 0.5, 1, 2, 4)) sweep$V <- vapply(sweep$c1, V_at, 0) sweep$ratio <- sweep$V / V_srs round(sweep, 4) c(break_even_c1 = break_even, as_fraction_of_c2 = break_even / c2) ``` ```{r scalars-3} #| echo: false be_frac <- break_even / c2 ``` The break-even point is a phase-one unit costing `r sprintf("%.0f", 100 * be_frac)` percent of a phase-two unit. Below that, two phases beat one; above it, you should have spent everything on the response. In this landscape the classification would have to be genuinely quick, a call from the vehicle rather than a walked transect. ```{r fig-cost} #| echo: false #| fig-cap: "Two-phase variance relative to a simple random sample, as the cost of a phase-one unit rises. The design pays for itself only to the left of the crossing point." #| fig-alt: "Rising curve of relative variance against phase-one cost, crossing the value one just before a phase-one cost of three." cs <- seq(0.1, 6, by = 0.05) cd2 <- data.frame(c1 = cs, ratio = vapply(cs, function(cc) V_at(cc) / V_srs, 0)) ggplot(cd2, aes(c1, ratio)) + geom_hline(yintercept = 1, colour = "#b5534e", linetype = "dashed", linewidth = 0.4) + geom_line(colour = "#275139", linewidth = 0.9) + annotate("point", x = break_even, y = 1, colour = "#b5534e", size = 3) + labs(title = "There is a price above which the first phase is a waste", subtitle = "Phase-two unit costs 10; dashed line is parity with a simple random sample", x = "cost of one phase-one unit", y = "variance relative to simple random sampling") + theme_te() ``` ## Two refinements worth knowing The subsampling fraction does not have to be the same in every stratum. Minimising the same variance with $\nu_h$ free gives the Neyman analogue, $\nu_h \propto S_h$, and the within term becomes $(\sum_h W_h S_h)^2/n_2$: ```{r optimal-nu} Sh <- sqrt(S2h) V_2ph_opt <- function(n1, n2) (1 - n1 / N) * S2 / n1 + (sum(W * Sh))^2 / n2 - within / n1 r2 <- sqrt((S2 - (sum(W * Sh))^2) * c2 / ((sum(W * Sh))^2 * c1)) n2b <- budget / (r2 * c1 + c2); n1b <- r2 * n2b c(n1 = n1b, n2 = n2b, V_uniform = V_2ph(n1_opt, n2_opt), V_optimal_nu = V_2ph_opt(n1b, n2b)) ``` ```{r scalars-4} #| echo: false nu_gain <- 1 - V_2ph_opt(n1b, n2b) / V_2ph(n1_opt, n2_opt) ``` Subsampling in proportion to the stratum standard deviations takes another `r sprintf("%.0f", 100 * nu_gain)` percent off the variance, for free, using information the first phase has already collected. It is the same Neyman argument one level down. The second refinement is a warning. The phase-one classification is allowed to be wrong. Because the estimator weights whatever classes the field call produced, misclassification does not bias anything; it only makes the classes separate the response less well, which shows up as a smaller between component. ```{r proxy-error} set.seed(3733) proxy <- do.call(rbind, lapply(c(0, 0.10, 0.20, 0.35), function(e) { cls <- stratum flip <- runif(N) < e cls[flip] <- sample(names(Nh), sum(flip), replace = TRUE) Wc <- as.numeric(table(factor(cls, levels = names(Nh)))) / N wi <- sum(Wc * tapply(y, cls, var)[names(Nh)]) r <- sqrt((S2 - wi) * c2 / (wi * c1)); n2e <- budget / (r * c1 + c2); n1e <- r * n2e Ve <- (1 - n1e / N) * S2 / n1e + wi * (1 / n2e - 1 / n1e) data.frame(error = e, within_share = wi / S2, V = Ve, ratio_to_srs = Ve / V_srs) })) round(proxy, 3) ``` ```{r scalars-5} #| echo: false r10 <- proxy$ratio_to_srs[proxy$error == 0.10] r35 <- proxy$ratio_to_srs[proxy$error == 0.35] ``` A tenth of the units misclassified moves the relative variance from `r sprintf("%.2f", proxy$ratio_to_srs[proxy$error == 0])` to `r sprintf("%.2f", r10)`, and at `r sprintf("%.0f", 100 * 0.35)` percent misclassification the design has slipped past parity to `r sprintf("%.2f", r35)`. Validity is safe, precision is not. A cheap classification that is cheap because it is careless will quietly spend your budget on nothing. ## The honest limit Everything here rests on phase one being a probability sample of the same frame as phase two, with phase two drawn inside it. Convenience units bolted on at either stage break the chain, and no amount of careful weighting repairs them. The variance formula is an approximation valid when the phase-one stratum counts are comfortably large. With a small first phase, or a class so rare that it collects a handful of units, the approximation runs optimistic in the wrong direction and the realised precision is worse than advertised. Two-phase designs are also more work to run than they look on paper: revisiting a specific subsample means keeping track of which units were which, in the field, months apart. ## Where to go next Three posts of design and no diagnostics yet. The last post in this group takes a finished stratified survey and asks the awkward questions: what happens when you stratify after the fact instead of before, what an out-of-date area figure does to the estimate, and how many strata are too many. ## References Neyman 1938 J Am Stat Assoc 33(201):101-116 (10.1080/01621459.1938.10503378) Rao 1973 Biometrika 60(1):125-133 (10.1093/biomet/60.1.125) Cochran 1977 Sampling Techniques, 3rd edn, Wiley; ISBN 978-0-471-16240-7 Sarndal, Swensson and Wretman 1992 Model Assisted Survey Sampling, Springer; ISBN 978-0-387-40620-6 Gregoire and Valentine 2008 Sampling Strategies for Natural Resources and the Environment, Chapman and Hall; ISBN 978-1-58488-370-8 ## Related tutorials - [Stratified random sampling in ecology](../stratified-random-sampling/) - [Allocating survey effort across strata](../allocating-effort-across-strata/) - [Checking a stratified design](../checking-a-stratified-design/) - [Unequal-probability spatial sampling](../unequal-probability-spatial-sampling/)