Stratified random sampling in ecology

Stratified random sampling in R from first principles: the weighted estimator, its design variance, and why splitting a landscape into strata sharpens a mean.

Walk any real landscape and the variation is not spread evenly. A wet hollow carries several times the biomass of the dry woodland floor next to it, and the two barely overlap. If you scatter plots at random over the whole area, most of the spread in your data comes from the fact that some plots landed in one habitat and some in the other. That spread is not news. You knew where the habitats were before you started.

Stratified sampling turns that prior knowledge into precision. You split the frame into strata, sample inside each one separately, and put the pieces back together with weights you already know. The estimator stays design based: no model, no assumption about the shape of anything, just the sampling plan you followed.

This post builds the estimator and its variance by hand, checks both against a simulation, and then shows the case where stratification buys nothing at all.

A frame you can count

The landscape is 1200 plots in three habitat strata. Their sizes are known from a habitat map, which is the whole point: the weights come from the frame, not from the sample.

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)
H       <- length(Nh)
stratum <- rep(names(Nh), Nh)
y       <- round(pmax(0, rnorm(N, rep(mu_h, Nh), rep(sd_h, Nh))), 2)

W   <- Nh / N                                    # stratum weights, known
mu  <- mean(y)                                   # population mean (we get to peek)
S2  <- var(y)                                    # population variance, divisor N - 1
S2h <- tapply(y, stratum, var)[names(Nh)]        # within-stratum variances
muh <- tapply(y, stratum, mean)[names(Nh)]
round(rbind(W = W, mean = muh, sd = sqrt(S2h)), 3)

     meadow wetland woodland
W     0.500   0.250    0.250
mean 11.985  26.079    5.138
sd    4.034   8.681    1.857

Figure 1: Plot-level biomass in the three strata. The wetland sits well above the meadow, the woodland well below, and the wetland is also the most variable stratum.

Two estimators for the same mean

A simple random sample of n plots estimates the mean with the plain sample mean. Its design variance is the textbook one, finite population correction included:

\[V_{\text{srs}}(\bar y) = \left(1 - \tfrac{n}{N}\right)\frac{S^2}{n}.\]

A stratified sample takes \(n_h\) plots inside stratum \(h\) and combines the stratum means with the known weights, \(\bar y_{\text{st}} = \sum_h W_h \bar y_h\). Because the strata are sampled independently, the variances simply add:

\[V_{\text{st}}(\bar y_{\text{st}}) = \sum_h W_h^2 \left(1 - \tfrac{n_h}{N_h}\right)\frac{S_h^2}{n_h}.\]

With a budget of 120 plots and effort split in proportion to stratum size:

n  <- 120
f  <- n / N
nh <- round(n * W)                                # proportional allocation
nh

  meadow  wetland woodland 
      60       30       30 

V_srs  <- (1 - f) * S2 / n
V_prop <- sum(W^2 * (1 - nh / Nh) * S2h / nh)
c(SE_srs = sqrt(V_srs), SE_stratified = sqrt(V_prop), ratio = V_srs / V_prop)

       SE_srs SE_stratified         ratio 
    0.8027738     0.4569169     3.0868246 

The stratified standard error is 0.457 against 0.803 for the simple random sample, a variance ratio of 3.09. To match the stratified design with random plots you would need about 371 plots instead of 120, for the same field season, on the same landscape, measuring the same thing. The only extra input was the habitat map.

Where the gain comes from

Split the population variance into a within part and a between part. The identity is exact for a finite population:

\[(N-1)S^2 = \sum_h (N_h - 1) S_h^2 + \sum_h N_h (\mu_h - \mu)^2.\]

within  <- sum((Nh - 1) * S2h) / (N - 1)
between <- sum(Nh * (muh - mu)^2) / (N - 1)
c(S2 = S2, within = within, between = between,
  check = within + between - S2, between_share = between / S2)

           S2        within       between         check between_share 
 8.592611e+01  2.778033e+01  5.814578e+01  8.526513e-14  6.766951e-01 

67.7 percent of the total variance lives between the strata. Now compare the two design variances under proportional allocation, where \(n_h/N_h\) equals \(n/N\) in every stratum:

\[V_{\text{srs}} - V_{\text{st}} = \frac{1-f}{n}\left(S^2 - \sum_h W_h S_h^2\right).\]

The bracket is the between-stratum component, up to a finite population bookkeeping term of order \(1/N\) (here the two versions differ by 0.056, against a total variance of 85.9). So the rule is short: stratification removes the between-stratum variance from your error, and nothing else. Everything else about the design follows from that one sentence.

Does the arithmetic survive a simulation?

Formulas for design variance are easy to write and easy to get wrong. Draw the samples and look.

idx_h <- split(seq_len(N), stratum)[names(Nh)]
B     <- 4000

draw_srs <- function() mean(y[sample.int(N, n)])
draw_str <- function() sum(W * vapply(seq_len(H), function(h)
                        mean(y[sample(idx_h[[h]], nh[h])]), 0))

set.seed(3711); mc_srs <- replicate(B, draw_srs())
set.seed(3712); mc_str <- replicate(B, draw_str())

rbind(simple     = c(mean = mean(mc_srs), variance = var(mc_srs), analytic = V_srs),
      stratified = c(mean = mean(mc_str), variance = var(mc_str), analytic = V_prop))

               mean  variance  analytic
simple     13.77937 0.6543158 0.6444459
stratified 13.79498 0.2070077 0.2087731

Both estimators centre on the true mean of 13.80, and both simulated variances land on their formulas: 0.654 against 0.644, and 0.207 against 0.209.

Figure 2: Sampling distributions of the two estimators over 4000 draws, same 120 plots of effort. Both are centred on the population mean; the stratified one is far narrower.

The variance estimator you actually report

In the field you do not know \(S_h^2\). You estimate it from the sample, stratum by stratum, and plug it in. That estimator is unbiased, and the interval it produces covers at the nominal rate:

one_survey <- function() {
  s   <- lapply(seq_len(H), function(h) y[sample(idx_h[[h]], nh[h])])
  est <- sum(W * vapply(s, mean, 0))
  v   <- sum(W^2 * (1 - nh / Nh) * vapply(s, var, 0) / nh)
  c(estimate = est, variance = v)
}

set.seed(3713)
surveys  <- t(replicate(B, one_survey()))
coverage <- mean(abs(surveys[, "estimate"] - mu) <= 1.96 * sqrt(surveys[, "variance"]))
c(mean_v_hat = mean(surveys[, "variance"]), true_V = V_prop, coverage = coverage)

mean_v_hat     true_V   coverage 
 0.2086173  0.2087731  0.9515000 

The average estimated variance is 0.2086 against a true 0.2088, and 95.2 percent of the intervals caught the mean. Nothing here leans on normality of the plot values; it leans on the central limit theorem applied inside each stratum, which is why very small \(n_h\) is a problem worth a post of its own.

When stratification buys nothing

The gain is the between-stratum variance. Take that away and the gain goes with it. Shift each stratum onto the common mean, leaving the within-stratum spread untouched, and re-run the same arithmetic:

y_flat <- y
for (h in names(Nh)) y_flat[stratum == h] <- y[stratum == h] - muh[h] + mu

S2_flat  <- var(y_flat)
S2h_flat <- tapply(y_flat, stratum, var)[names(Nh)]
c(V_srs = (1 - f) * S2_flat / n,
  V_str = sum(W^2 * (1 - nh / Nh) * S2h_flat / nh),
  ratio = ((1 - f) * S2_flat / n) / sum(W^2 * (1 - nh / Nh) * S2h_flat / nh))

    V_srs     V_str     ratio 
0.2083525 0.2087731 0.9979854 

The ratio is 0.998. The wetland is still 4.7 times as variable as the woodland, the map is still correct, the allocation is still proportional, and the design has gained nothing measurable. Strata that do not separate the response are decoration.

Sweeping the separation from zero up to one and a half times the real gaps shows the whole curve:

k_grid <- seq(0, 1.5, by = 0.125)
ratio_k <- vapply(k_grid, function(k) {
  yk <- y
  for (h in names(Nh))
    yk[stratum == h] <- y[stratum == h] - muh[h] + mu + k * (muh[h] - mu)
  S2k  <- var(yk)
  S2hk <- tapply(yk, stratum, var)[names(Nh)]
  ((1 - f) * S2k / n) / sum(W^2 * (1 - nh / Nh) * S2hk / nh)
}, 0)
setNames(round(ratio_k, 2), k_grid)

    0 0.125  0.25 0.375   0.5 0.625  0.75 0.875     1 1.125  1.25 1.375   1.5 
 1.00  1.03  1.13  1.29  1.52  1.81  2.17  2.60  3.09  3.64  4.26  4.95  5.70 

Figure 3: Variance ratio against how far apart the stratum means are, as a multiple of the real separation. At zero separation the ratio is one: identical within-stratum spreads, identical precision.

The honest limit

Stratification is optimal for the variable you stratified on, and for no other variable automatically. A partition that separates biomass may cut straight across the gradient that drives species richness, in which case the richness estimate gains nothing while the biomass estimate gains a factor of 3.1. Surveys usually carry more than one response, and the strata are one decision serving all of them.

Two further conditions do real work and both are about the frame rather than the statistics. The stratum sizes \(N_h\) have to be the true sizes of the sets you sampled from, because they enter the estimator directly. And every plot has to belong to exactly one stratum with a known, non-zero chance of selection. Get either wrong and the estimator stops being unbiased, which is a different kind of problem from being imprecise.

Where to go next

Proportional allocation is the safe default, not the best one. The wetland here is 4.7 times as variable as the woodland, and sending equal fractions of effort to both is leaving precision on the table. The next post works out the allocation that minimises the variance, what it costs when a stratum is expensive to reach, and how much of the promised gain survives being planned from a small pilot.

References

Neyman 1934 J R Stat Soc 97(4):558-625 (10.2307/2342192)

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

Thompson 2012 Sampling, 3rd edn, Wiley; ISBN 978-0-470-40231-3

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

Related tutorials

• Allocating survey effort across strata
• Two-phase sampling when strata are unknown
• Checking a stratified design
• Spatially balanced sampling with GRTS