Checking a survey design

survey design

model diagnostics

monitoring

ecology tutorial

Three diagnostics for a spatially balanced ecological survey in R: realised inclusion probabilities, systematic grid aliasing, and what the frame leaves out.

Author

Tidy Ecology

Published

2026-07-25

The three previous posts built a spatially balanced design, measured what it buys, and showed how to weight when the inclusion probabilities are unequal. This one asks what can still go wrong. Survey design is unusual among the topics on this blog in that the assumptions are mostly under your control: you chose the frame, you set the inclusion probabilities, you decided how the sites are laid out. That makes the failures administrative rather than statistical, and administrative failures are quiet. Three checks are worth running before a monitoring programme goes to the field, and the third one cannot be run at all after the fact.

library(ggplot2); library(MASS)
theme_te <- function() {
  theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e6da", linewidth = 0.3),
          plot.title = element_text(colour = "#16241d", face = "bold", size = 13),
          plot.subtitle = element_text(colour = "#5d6b61", size = 10),
          axis.title = element_text(colour = "#46604a"),
          axis.text = element_text(colour = "#5d6b61"),
          legend.title = element_text(colour = "#46604a"),
          legend.text = element_text(colour = "#5d6b61"))
}
pal_forest <- "#275139"; pal_gold <- "#cda23f"; pal_rust <- "#b5534e"

grts_address <- function(x, y, L = 8L) {
  n <- length(x); addr <- numeric(n)
  recurse <- function(idx, xlo, xhi, ylo, yhi, level, prefix) {
    if (level >= L || length(idx) <= 1L) { addr[idx] <<- prefix; return(invisible(NULL)) }
    xm <- (xlo + xhi)/2; ym <- (ylo + yhi)/2
    q <- (x[idx] >= xm) + 2L*(y[idx] >= ym); perm <- sample(0:3)
    for (qq in 0:3) {
      sub <- idx[q == qq]; if (length(sub) == 0L) next
      cxlo <- if (qq %% 2L == 1L) xm else xlo; cxhi <- if (qq %% 2L == 1L) xhi else xm
      cylo <- if (qq >= 2L) ym else ylo;       cyhi <- if (qq >= 2L) yhi else ym
      recurse(sub, cxlo, cxhi, cylo, cyhi, level + 1L, prefix*4 + perm[qq + 1L])
    }
  }
  recurse(seq_len(n), 0, 1, 0, 1, 0L, 0); addr
}
grts_draw <- function(x, y, m, pip = NULL, L = 8L) {
  n <- length(x); if (is.null(pip)) pip <- rep(m/n, n)
  a <- grts_address(x, y, L); ord <- order(a, runif(n))
  cs <- cumsum(pip[ord]); picks <- runif(1) + 0:(m - 1)
  ord[findInterval(picks, c(0, cs), rightmost.closed = TRUE)]
}

Check one: do the realised inclusion probabilities match the plan?

Everything downstream rests on pi_i. The Horvitz-Thompson weights are the reciprocals of the numbers you wrote down, so if the draw does not actually select units at those rates, the weighting corrects for a design that never happened. This is worth verifying because it is easy to break: capping probabilities at one, renormalising after the cap, dropping units from the frame, or a bug in the size measure all shift the realised rates away from the intended ones.

The check is a simulation of your own design. Run the draw many times against the frame you plan to use, count how often each unit is selected, and compare with the intended probabilities.

set.seed(404)
side <- 22
g <- expand.grid(x = seq(0.02, 0.98, length.out = side),
                 y = seq(0.02, 0.98, length.out = side))
N <- nrow(g); m <- 40
w <- 1 + 3 * g$x                                   # intended density rises to the east
pip <- m * w / sum(w); pip <- pmin(pip, 0.999); pip <- m * pip / sum(pip)

set.seed(111)
R <- 6000
counts <- integer(N)
for (r in seq_len(R)) { s <- grts_draw(g$x, g$y, m, pip = pip); counts[s] <- counts[s] + 1L }
pi_emp <- counts / R
cor_pi <- cor(pi_emp, pip); dev_pi <- max(abs(pi_emp - pip))

Over 6000 draws the realised and intended inclusion probabilities correlate at 0.9926, with a largest absolute deviation of 0.0116, comfortably inside Monte Carlo error for this many replicates. The design does what the specification says. A failed version of this check looks obvious: a systematic gap between the two, or a set of units pinned at zero because they fell out of the frame.

A scatter plot of realised against intended inclusion probability. Points lie tightly along a dashed one-to-one line rising from lower left to upper right. — Figure 1: Realised against intended inclusion probabilities over six thousand draws. Points sit on the identity line, so the implemented design matches its specification.

Check two: would a systematic grid have aliased?

A regular grid is the other way to guarantee even spatial coverage, and it is tempting because it is simple to explain and simple to navigate to. Its weakness is that it has almost no randomisation left: a grid with a random start has only as many distinct configurations as there are offsets. If the landscape carries periodic structure at the grid spacing, every site lands at the same phase of the cycle, and the sample tells you about that phase and nothing else. The estimate stays roughly unbiased across random starts, but its variance explodes.

Periodic structure in ecology is not exotic. Agricultural furrows, planted rows, dune and swale systems, ridge and slack, drainage ditches and terraces all impose regular spacing on a landscape.

sideF <- 30
gf <- expand.grid(r = 1:sideF, c = 1:sideF)
gf$x <- (gf$c - 0.5)/sideF; gf$y <- (gf$r - 0.5)/sideF
Nf <- nrow(gf); step <- sideF/6; base_rows <- seq(1, sideF, by = step)

set.seed(707)
per <- 6                                          # waves match the grid spacing
y_per <- 10 + 3*sin(2*pi*per*gf$x) + rnorm(Nf, 0, 0.5)
Dmf <- as.matrix(dist(cbind(gf$x, gf$y)))
set.seed(909)
y_smooth <- as.numeric(mvrnorm(1, rep(10, Nf), 3*exp(-Dmf/0.2)))

sys_sample <- function(off) {
  rr <- ((base_rows - 1 + off[1]) %% sideF) + 1
  cc <- ((base_rows - 1 + off[2]) %% sideF) + 1
  which(gf$r %in% rr & gf$c %in% cc)
}

set.seed(808)
Rr <- 3000
sys_per <- numeric(Rr); grts_per <- numeric(Rr)
sys_sm  <- numeric(Rr); grts_sm  <- numeric(Rr)
for (k in seq_len(Rr)) {
  ss <- sys_sample(sample(0:(step - 1), 2, replace = TRUE))
  sg <- grts_draw(gf$x, gf$y, length(ss))
  sys_per[k] <- mean(y_per[ss]);    grts_per[k] <- mean(y_per[sg])
  sys_sm[k]  <- mean(y_smooth[ss]); grts_sm[k]  <- mean(y_smooth[sg])
}
ratio_per <- var(sys_per) / var(grts_per)
ratio_sm  <- var(sys_sm)  / var(grts_sm)

On the periodic landscape the systematic grid has 40 times the variance of GRTS at the same sample size: standard deviation 2.153 against 0.342. Any single systematic survey returns a confident, precise-looking number that is a property of the phase it happened to land on. On the smooth landscape, where no periodicity exists to lock onto, the two designs are comparable and GRTS is modestly ahead, with a variance ratio of 0.92.

The asymmetry is the point. GRTS is never much worse than a grid and is enormously better in the case that ruins a grid, and you do not need to know in advance whether your landscape has periodic structure at your spacing. That is what the randomised hierarchy is protecting you from.

Two density curves of the estimated mean. The GRTS curve is a narrow spike near the true mean; the systematic-grid curve is broad and flat, spanning several units either side. — Figure 2: Sampling distributions of the estimated mean on a periodic landscape whose spacing matches the grid. The systematic design spreads wildly across random starts; GRTS is unaffected.

Check three: what did the frame leave out?

The first two checks can be run from your desk. This one cannot be run from the data at all, and it is the one that most often decides whether a monitoring programme means what its report says.

A design-based estimate is unbiased for the frame. If places were quietly dropped because they are steep, private, far from a road or under water at the time of the survey, then the frame is a biased subset of the population, and every standard error you compute is honest about sampling error while silent about the omission.

set.seed(707)
Dm <- as.matrix(dist(g))
elev <- 0.6*g$x + 0.8*g$y                          # ground rises to the north-east
z <- 8 + 5*elev + as.numeric(mvrnorm(1, rep(0, N), 1.5*exp(-Dm/0.2)))
access <- elev < quantile(elev, 0.82)              # the steep corner is never surveyed
mu_pop <- mean(z); mu_frame <- mean(z[access])

set.seed(808)
idx <- which(access)
est <- replicate(3000, mean(z[idx[grts_draw(g$x[idx], g$y[idx], m)]]))
err_frame <- mean(est) - mu_frame
err_pop   <- mean(est) - mu_pop
se_design <- sd(est)
n_se <- abs(err_pop) / se_design

The frame covers 82% of the landscape, leaving out the steepest fifth. Against the frame mean the design is exactly as advertised: average error +0.000. Against the population mean it is off by -0.562, which is 3.3 design standard errors. No diagnostic run on the sample can see this, because every sampled unit is a legitimate member of the frame and the estimator is behaving correctly. The error lives in the gap between the list you sampled and the population you are going to write about.

Left panel: a square landscape with the north-east corner marked as excluded from the frame. Right panel: a density curve centred on a dashed frame-mean line, with a separate solid line marking the higher population mean well outside the curve. — Figure 3: The omitted region and its consequence. Left: the frame excludes the steep north-east corner. Right: the sampling distribution is centred on the frame mean, not on the population mean.

What the checks can and cannot settle

The first two checks are genuine tests with a possible failure: the realised probabilities can disagree with the plan, and a grid can alias. Run them before the field season, when the answer can still change the design.

The third is not a test. It is a statement you have to write down and defend, because the data will never contradict it. The useful discipline is to record what the frame excludes and why, and to say plainly whether the reported quantity is a property of the surveyed population or of the accessible part of it. A design-based survey gives you an unusually clean guarantee, and the guarantee is precisely as wide as the frame.

References

Stevens, D.L. and Olsen, A.R. (2004). Journal of the American Statistical Association 99(465):262-278 (10.1198/016214504000000250).

Stevens, D.L. and Olsen, A.R. (2003). Environmetrics 14(6):593-610 (10.1002/env.606).

Olsen, A.R., Sedransk, J., Edwards, D., Gotway, C.A., Liggett, W., Rathbun, S., Reckhow, K.H. and Young, L.J. (1999). Environmental Monitoring and Assessment 54(1):1-45 (10.1023/A:1005823911258).

Theobald, D.M., Stevens, D.L., White, D., Urquhart, N.S., Olsen, A.R. and Norman, J.B. (2007). Environmental Management 40(1):134-146 (10.1007/s00267-005-0199-x).

--- title: "Checking a survey design" description: "Three diagnostics for a spatially balanced ecological survey in R: realised inclusion probabilities, systematic grid aliasing, and what the frame leaves out." date: "2026-07-25 12:00" categories: [R, survey design, model diagnostics, monitoring, ecology tutorial] image: thumbnail.png image-alt: "Three small diagnostic panels: realised against intended inclusion probabilities, sampling distributions under a systematic grid and GRTS, and a frame with an omitted corner." --- The three previous posts built a spatially balanced design, measured what it buys, and showed how to weight when the inclusion probabilities are unequal. This one asks what can still go wrong. Survey design is unusual among the topics on this blog in that the assumptions are mostly under your control: you chose the frame, you set the inclusion probabilities, you decided how the sites are laid out. That makes the failures administrative rather than statistical, and administrative failures are quiet. Three checks are worth running before a monitoring programme goes to the field, and the third one cannot be run at all after the fact. ```{r setup} #| message: false #| warning: false library(ggplot2); library(MASS) theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e6da", linewidth = 0.3), plot.title = element_text(colour = "#16241d", face = "bold", size = 13), plot.subtitle = element_text(colour = "#5d6b61", size = 10), axis.title = element_text(colour = "#46604a"), axis.text = element_text(colour = "#5d6b61"), legend.title = element_text(colour = "#46604a"), legend.text = element_text(colour = "#5d6b61")) } pal_forest <- "#275139"; pal_gold <- "#cda23f"; pal_rust <- "#b5534e" grts_address <- function(x, y, L = 8L) { n <- length(x); addr <- numeric(n) recurse <- function(idx, xlo, xhi, ylo, yhi, level, prefix) { if (level >= L || length(idx) <= 1L) { addr[idx] <<- prefix; return(invisible(NULL)) } xm <- (xlo + xhi)/2; ym <- (ylo + yhi)/2 q <- (x[idx] >= xm) + 2L*(y[idx] >= ym); perm <- sample(0:3) for (qq in 0:3) { sub <- idx[q == qq]; if (length(sub) == 0L) next cxlo <- if (qq %% 2L == 1L) xm else xlo; cxhi <- if (qq %% 2L == 1L) xhi else xm cylo <- if (qq >= 2L) ym else ylo; cyhi <- if (qq >= 2L) yhi else ym recurse(sub, cxlo, cxhi, cylo, cyhi, level + 1L, prefix*4 + perm[qq + 1L]) } } recurse(seq_len(n), 0, 1, 0, 1, 0L, 0); addr } grts_draw <- function(x, y, m, pip = NULL, L = 8L) { n <- length(x); if (is.null(pip)) pip <- rep(m/n, n) a <- grts_address(x, y, L); ord <- order(a, runif(n)) cs <- cumsum(pip[ord]); picks <- runif(1) + 0:(m - 1) ord[findInterval(picks, c(0, cs), rightmost.closed = TRUE)] } ``` ## Check one: do the realised inclusion probabilities match the plan? Everything downstream rests on `pi_i`. The Horvitz-Thompson weights are the reciprocals of the numbers you wrote down, so if the draw does not actually select units at those rates, the weighting corrects for a design that never happened. This is worth verifying because it is easy to break: capping probabilities at one, renormalising after the cap, dropping units from the frame, or a bug in the size measure all shift the realised rates away from the intended ones. The check is a simulation of your own design. Run the draw many times against the frame you plan to use, count how often each unit is selected, and compare with the intended probabilities. ```{r check-pi} set.seed(404) side <- 22 g <- expand.grid(x = seq(0.02, 0.98, length.out = side), y = seq(0.02, 0.98, length.out = side)) N <- nrow(g); m <- 40 w <- 1 + 3 * g$x # intended density rises to the east pip <- m * w / sum(w); pip <- pmin(pip, 0.999); pip <- m * pip / sum(pip) set.seed(111) R <- 6000 counts <- integer(N) for (r in seq_len(R)) { s <- grts_draw(g$x, g$y, m, pip = pip); counts[s] <- counts[s] + 1L } pi_emp <- counts / R cor_pi <- cor(pi_emp, pip); dev_pi <- max(abs(pi_emp - pip)) ``` Over `r R` draws the realised and intended inclusion probabilities correlate at `r sprintf("%.4f", cor_pi)`, with a largest absolute deviation of `r sprintf("%.4f", dev_pi)`, comfortably inside Monte Carlo error for this many replicates. The design does what the specification says. A failed version of this check looks obvious: a systematic gap between the two, or a set of units pinned at zero because they fell out of the frame. ```{r fig-pi} #| echo: false #| fig-cap: "Realised against intended inclusion probabilities over six thousand draws. Points sit on the identity line, so the implemented design matches its specification." #| fig-alt: "A scatter plot of realised against intended inclusion probability. Points lie tightly along a dashed one-to-one line rising from lower left to upper right." ggplot(data.frame(intended = pip, realised = pi_emp), aes(intended, realised)) + geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "#5d6b61") + geom_point(colour = pal_forest, alpha = 0.5, size = 1.6) + labs(title = "The draw selects units at the rates you specified", subtitle = "Dashed line: exact agreement.", x = "Intended inclusion probability", y = "Realised inclusion probability") + theme_te() ``` ## Check two: would a systematic grid have aliased? A regular grid is the other way to guarantee even spatial coverage, and it is tempting because it is simple to explain and simple to navigate to. Its weakness is that it has almost no randomisation left: a grid with a random start has only as many distinct configurations as there are offsets. If the landscape carries periodic structure at the grid spacing, every site lands at the same phase of the cycle, and the sample tells you about that phase and nothing else. The estimate stays roughly unbiased across random starts, but its variance explodes. Periodic structure in ecology is not exotic. Agricultural furrows, planted rows, dune and swale systems, ridge and slack, drainage ditches and terraces all impose regular spacing on a landscape. ```{r check-alias} sideF <- 30 gf <- expand.grid(r = 1:sideF, c = 1:sideF) gf$x <- (gf$c - 0.5)/sideF; gf$y <- (gf$r - 0.5)/sideF Nf <- nrow(gf); step <- sideF/6; base_rows <- seq(1, sideF, by = step) set.seed(707) per <- 6 # waves match the grid spacing y_per <- 10 + 3*sin(2*pi*per*gf$x) + rnorm(Nf, 0, 0.5) Dmf <- as.matrix(dist(cbind(gf$x, gf$y))) set.seed(909) y_smooth <- as.numeric(mvrnorm(1, rep(10, Nf), 3*exp(-Dmf/0.2))) sys_sample <- function(off) { rr <- ((base_rows - 1 + off[1]) %% sideF) + 1 cc <- ((base_rows - 1 + off[2]) %% sideF) + 1 which(gf$r %in% rr & gf$c %in% cc) } set.seed(808) Rr <- 3000 sys_per <- numeric(Rr); grts_per <- numeric(Rr) sys_sm <- numeric(Rr); grts_sm <- numeric(Rr) for (k in seq_len(Rr)) { ss <- sys_sample(sample(0:(step - 1), 2, replace = TRUE)) sg <- grts_draw(gf$x, gf$y, length(ss)) sys_per[k] <- mean(y_per[ss]); grts_per[k] <- mean(y_per[sg]) sys_sm[k] <- mean(y_smooth[ss]); grts_sm[k] <- mean(y_smooth[sg]) } ratio_per <- var(sys_per) / var(grts_per) ratio_sm <- var(sys_sm) / var(grts_sm) ``` On the periodic landscape the systematic grid has `r sprintf("%.0f", ratio_per)` times the variance of GRTS at the same sample size: standard deviation `r sprintf("%.3f", sd(sys_per))` against `r sprintf("%.3f", sd(grts_per))`. Any single systematic survey returns a confident, precise-looking number that is a property of the phase it happened to land on. On the smooth landscape, where no periodicity exists to lock onto, the two designs are comparable and GRTS is modestly ahead, with a variance ratio of `r sprintf("%.2f", ratio_sm)`. The asymmetry is the point. GRTS is never much worse than a grid and is enormously better in the case that ruins a grid, and you do not need to know in advance whether your landscape has periodic structure at your spacing. That is what the randomised hierarchy is protecting you from. ```{r fig-alias} #| echo: false #| fig-cap: "Sampling distributions of the estimated mean on a periodic landscape whose spacing matches the grid. The systematic design spreads wildly across random starts; GRTS is unaffected." #| fig-alt: "Two density curves of the estimated mean. The GRTS curve is a narrow spike near the true mean; the systematic-grid curve is broad and flat, spanning several units either side." da <- rbind(data.frame(est = sys_per, design = "Systematic grid"), data.frame(est = grts_per, design = "GRTS")) ggplot(da, aes(est, fill = design, colour = design)) + geom_density(alpha = 0.35, linewidth = 0.7) + geom_vline(xintercept = mean(y_per), linetype = "dashed", colour = "#5d6b61") + scale_fill_manual(values = c("GRTS" = pal_forest, "Systematic grid" = pal_rust)) + scale_colour_manual(values = c("GRTS" = pal_forest, "Systematic grid" = pal_rust)) + labs(title = "A grid that matches the landscape period fails badly", subtitle = "Dashed line: the true frame mean.", x = "Estimated mean", y = "Density", fill = NULL, colour = NULL) + theme_te() ``` ## Check three: what did the frame leave out? The first two checks can be run from your desk. This one cannot be run from the data at all, and it is the one that most often decides whether a monitoring programme means what its report says. A design-based estimate is unbiased **for the frame**. If places were quietly dropped because they are steep, private, far from a road or under water at the time of the survey, then the frame is a biased subset of the population, and every standard error you compute is honest about sampling error while silent about the omission. ```{r check-frame} set.seed(707) Dm <- as.matrix(dist(g)) elev <- 0.6*g$x + 0.8*g$y # ground rises to the north-east z <- 8 + 5*elev + as.numeric(mvrnorm(1, rep(0, N), 1.5*exp(-Dm/0.2))) access <- elev < quantile(elev, 0.82) # the steep corner is never surveyed mu_pop <- mean(z); mu_frame <- mean(z[access]) set.seed(808) idx <- which(access) est <- replicate(3000, mean(z[idx[grts_draw(g$x[idx], g$y[idx], m)]])) err_frame <- mean(est) - mu_frame err_pop <- mean(est) - mu_pop se_design <- sd(est) n_se <- abs(err_pop) / se_design ``` The frame covers `r sprintf("%.0f%%", 100*mean(access))` of the landscape, leaving out the steepest fifth. Against the frame mean the design is exactly as advertised: average error `r sprintf("%+.3f", err_frame)`. Against the population mean it is off by `r sprintf("%+.3f", err_pop)`, which is `r sprintf("%.1f", n_se)` design standard errors. No diagnostic run on the sample can see this, because every sampled unit is a legitimate member of the frame and the estimator is behaving correctly. The error lives in the gap between the list you sampled and the population you are going to write about. ```{r fig-frame} #| echo: false #| fig-cap: "The omitted region and its consequence. Left: the frame excludes the steep north-east corner. Right: the sampling distribution is centred on the frame mean, not on the population mean." #| fig-alt: "Left panel: a square landscape with the north-east corner marked as excluded from the frame. Right panel: a density curve centred on a dashed frame-mean line, with a separate solid line marking the higher population mean well outside the curve." p1 <- ggplot(data.frame(g, access = access), aes(x, y, colour = access)) + geom_point(size = 1.8) + scale_colour_manual(values = c("TRUE" = pal_forest, "FALSE" = pal_rust), labels = c("omitted", "in frame"), name = NULL) + coord_equal() + labs(title = "The frame", x = NULL, y = NULL) + theme_te() + theme(axis.text = element_blank(), legend.position = "bottom") p2 <- ggplot(data.frame(est = est), aes(est)) + geom_density(fill = pal_forest, colour = pal_forest, alpha = 0.35, linewidth = 0.7) + geom_vline(xintercept = mu_frame, linetype = "dashed", colour = "#5d6b61") + geom_vline(xintercept = mu_pop, colour = pal_rust, linewidth = 0.8) + labs(title = "The estimate", subtitle = "Dashed: frame mean. Red: population mean.", x = "Estimated mean", y = "Density") + coord_cartesian(xlim = range(c(est, mu_pop)) + c(-0.1, 0.1)) + theme_te() grid::grid.newpage() grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(p1, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(p2, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## What the checks can and cannot settle The first two checks are genuine tests with a possible failure: the realised probabilities can disagree with the plan, and a grid can alias. Run them before the field season, when the answer can still change the design. The third is not a test. It is a statement you have to write down and defend, because the data will never contradict it. The useful discipline is to record what the frame excludes and why, and to say plainly whether the reported quantity is a property of the surveyed population or of the accessible part of it. A design-based survey gives you an unusually clean guarantee, and the guarantee is precisely as wide as the frame. ## Related tutorials - [Spatially balanced sampling with GRTS](../spatially-balanced-sampling-and-grts/) - [Design-based variance and spatial balance](../design-based-variance-and-spatial-balance/) - [Unequal-probability spatial sampling](../unequal-probability-spatial-sampling/) - [Survey design and detection visits](../survey-design-detection-visits/) ## References Stevens, D.L. and Olsen, A.R. (2004). Journal of the American Statistical Association 99(465):262-278 (10.1198/016214504000000250). Stevens, D.L. and Olsen, A.R. (2003). Environmetrics 14(6):593-610 (10.1002/env.606). Olsen, A.R., Sedransk, J., Edwards, D., Gotway, C.A., Liggett, W., Rathbun, S., Reckhow, K.H. and Young, L.J. (1999). Environmental Monitoring and Assessment 54(1):1-45 (10.1023/A:1005823911258). Theobald, D.M., Stevens, D.L., White, D., Urquhart, N.S., Olsen, A.R. and Norman, J.B. (2007). Environmental Management 40(1):134-146 (10.1007/s00267-005-0199-x).