Checking a monitoring design

Three checks on a monitoring plan in R: what a short pilot does to power, which kind of site dropout biases a trend, and why observer drift is invisible.

The previous three posts built a monitoring design and gave it a number: this programme has such and such power to detect such and such a decline. That number rests on assumptions, and this post attacks three of them with simulation.

The first is that the variance components are known. The second is that the sites measured in year twenty are the sites measured in year one. The third is that a count means the same thing in both years. Each check produces a number, and none of them is reassuring.

The setup

The planned programme: 30 permanent sites, two counts per site per year, twenty years, testing for a three per cent annual decline. The power calculation from the previous post says this design works.

sd_site <- 0.55; sd_year <- 0.20; sd_sy <- 0.30; sd_res <- 0.35
beta_true <- -0.03
n_rep <- 2

se_trend <- function(n_year, n_site, v = c(sd_year^2, sd_sy^2, sd_res^2), k = n_rep) {
  Sxx <- n_year * (n_year^2 - 1) / 12
  sqrt((v[1] + (v[2] + v[3] / k) / n_site) / Sxx)
}
trend_power <- function(n_year, n_site, beta, v = c(sd_year^2, sd_sy^2, sd_res^2), alpha = 0.05) {
  df <- n_year - 2; ncp <- abs(beta) / se_trend(n_year, n_site, v)
  crit <- qt(1 - alpha / 2, df)
  pt(-crit, df, ncp) + (1 - pt(crit, df, ncp))
}
true_power <- trend_power(20, 30, beta_true)
round(c(se = se_trend(20, 30), power = true_power), 4)

    se  power 
0.0082 0.9312 

sim_monitoring <- function(n_site, n_year, n_rep, beta) {
  s  <- rnorm(n_site, 0, sd_site)
  a  <- rnorm(n_year, 0, sd_year)
  tc <- (1:n_year) - mean(1:n_year)
  sa <- matrix(rnorm(n_site * n_year, 0, sd_sy), n_site, n_year)
  d  <- expand.grid(rep = seq_len(n_rep), site = seq_len(n_site), year = seq_len(n_year))
  d$y <- 3 + beta * tc[d$year] + s[d$site] + a[d$year] +
         sa[cbind(d$site, d$year)] + rnorm(nrow(d), 0, sd_res)
  d$site <- factor(d$site); d$yearf <- factor(d$year); d$t <- d$year
  d
}
vc_anova <- function(d, n_site, n_year, n_rep) {
  y <- d$y; gm <- mean(y)
  m_si <- tapply(y, d$site, mean); m_yr <- tapply(y, d$yearf, mean)
  m_cell <- tapply(y, list(d$site, d$yearf), mean)
  SS_year <- n_site * n_rep * sum((m_yr - gm)^2)
  SS_sy   <- n_rep * sum((m_cell - outer(m_si, m_yr, "+") + gm)^2)
  SS_res  <- sum((y - m_cell[cbind(as.integer(d$site), as.integer(d$yearf))])^2)
  MS_year <- SS_year / (n_year - 2)                  # one year df spent on the trend
  MS_sy   <- SS_sy / ((n_site - 1) * (n_year - 1))
  MS_res  <- SS_res / (n_site * n_year * (n_rep - 1))
  c(year = (MS_year - MS_sy) / (n_rep * n_site),
    sy   = (MS_sy - MS_res) / n_rep,
    res  = MS_res)
}

Check one: the power number is an estimate

The variance components come from somewhere, and that somewhere is usually a pilot: a couple of seasons at a handful of sites, run before the real programme starts. A pilot estimates the year component from a handful of years, and the year component is the one that dominates trend precision.

Here are 2,000 pilots of five years at ten sites, each feeding its own components into the same power calculation for the same planned design.

set.seed(370)
pilot_power <- t(replicate(2000, {
  d <- sim_monitoring(10, 5, n_rep, beta_true)
  d$y <- residuals(lm(y ~ t, data = d))
  v <- vc_anova(d, 10, 5, n_rep)
  c(v_year = unname(v["year"]),
    power  = trend_power(20, 30, beta_true, pmax(unname(v), 0)))
}))
round(c(true_year_component = sd_year^2,
        median_estimate     = median(pilot_power[, "v_year"]),
        negative_estimates  = mean(pilot_power[, "v_year"] < 0),
        true_power          = true_power,
        median_power        = median(pilot_power[, "power"]),
        lower_5pc           = unname(quantile(pilot_power[, "power"], 0.05)),
        pc_claiming_over_95 = mean(pilot_power[, "power"] > 0.95)), 4)

true_year_component     median_estimate  negative_estimates          true_power 
             0.0400              0.0293              0.1595              0.9312 
       median_power           lower_5pc pc_claiming_over_95 
             0.9768              0.5229              0.5625 

The true power is 0.931. The median pilot claims 0.977, and 56 per cent of pilots claim above 0.95. One pilot in twenty produces a plan whose real power is fine but whose own calculation says 0.52 or less.

The bias has a mechanism. Five years give roughly three degrees of freedom for the year component, its moment estimator goes negative in 16 per cent of pilots and is truncated at zero, and its median is 0.0293 against a true 0.04. A year component estimated too small makes the design look better than it is, and it is estimated too small more often than not.

Figure 1: Power claimed by 2,000 five-year pilots for the same twenty-year design. The dashed line is the true power of that design. The distribution is not centred on it and it is not narrow: the mass piled against one comes from pilots that estimated the year component at or below zero.

The practical response is not to abandon the pilot but to report the calculation as an interval, by running the power formula over the plausible range of the year component rather than at its point estimate. A design that only works at the optimistic end is not a design.

Check two: which kind of dropout matters

Sites are lost. Access is withdrawn, a landowner changes, a volunteer stops. The question is not whether attrition biases the trend but which kind does.

run_attrition <- function(mode, B = 1500, n_year = 15, n_site = 30, rate = 0.08) {
  tc <- (1:n_year) - mean(1:n_year)
  t(replicate(B, {
    s <- rnorm(n_site, 0, sd_site); a <- rnorm(n_year, 0, sd_year)
    alive <- rep(TRUE, n_site); annual <- numeric(n_year)
    for (j in 1:n_year) {
      if (j > 1 && mode != "none") {
        p_drop <- if (mode == "random") rep(rate, n_site)
                  else 2 * rate / (1 + exp(3 * (s - mean(s))))   # poor sites go first
        alive <- alive & !(runif(n_site) < p_drop)
      }
      ids <- which(alive)
      annual[j] <- mean(3 + beta_true * tc[j] + s[ids] + a[j] +
                        rnorm(length(ids), 0, sd_sy) +
                        rnorm(length(ids), 0, sd_res / sqrt(n_rep)))
    }
    c(slope = unname(coef(lm(annual ~ tc))[2]), left = sum(alive))
  }))
}
set.seed(3701)
att <- t(sapply(c("none", "random", "selective"), function(m) {
  r <- run_attrition(m)
  c(mean_slope = mean(r[, "slope"]), bias = mean(r[, "slope"]) - beta_true,
    sites_left = mean(r[, "left"]))
}))
round(att, 4)

          mean_slope   bias sites_left
none         -0.0298 0.0002    30.0000
random       -0.0294 0.0006     9.3327
selective    -0.0054 0.0246    11.4380

Random attrition removes two thirds of the network over fifteen years, leaving 9.3 of the original 30 sites, and biases the trend by +0.0006 against a true -0.03. It costs precision and nothing else.

Attrition that depends on the site itself is a different matter. When the sites most likely to be lost are the poor ones, the surviving network improves for reasons that have nothing to do with the population. The same overall loss rate leaves the estimated trend at -0.0054, which is 82 per cent of the true decline erased.

Figure 2: Estimated trend under three attrition regimes, all losing sites at the same rate. Random loss sits on the truth. Loss that depends on the state of the site moves the estimate most of the way to zero, and nothing in the surviving data marks it as a problem.

The diagnostic is available before the analysis: compare the early values of the sites that survived with the early values of the sites that were lost. If those two groups differ, the surviving network is not the network the programme started with, and the trend is partly a trend in membership.

Check three: the thing the design cannot see

A count is a population multiplied by a detection probability. A monitoring programme measures the product, and the trend it reports is the trend in the product.

drift_run <- function(pop_trend, detection_trend, B = 1500, n_year = 20, n_site = 30) {
  tc <- (1:n_year) - mean(1:n_year)
  t(replicate(B, {
    d <- sim_monitoring(n_site, n_year, n_rep, pop_trend)
    d$y <- d$y + detection_trend * tc[d$year]      # observers slowly change
    annual <- tapply(d$y, d$yearf, mean)
    co <- summary(lm(annual ~ tc))$coefficients[2, ]
    c(est = unname(co[1]), sig = as.numeric(abs(co[1] / co[2]) > qt(0.975, n_year - 2)))
  }))
}
set.seed(3702)
flat <- drift_run(0, -0.015)
decl <- drift_run(-0.03, -0.015)
round(c(stable_pop_mean_estimate = mean(flat[, "est"]),
        stable_pop_pc_significant = 100 * mean(flat[, "sig"]),
        declining_pop_mean_estimate = mean(decl[, "est"]),
        sum_of_the_two_trends = -0.03 - 0.015), 4)

   stable_pop_mean_estimate   stable_pop_pc_significant 
                    -0.0147                     39.2667 
declining_pop_mean_estimate       sum_of_the_two_trends 
                    -0.0449                     -0.0450 

A population that is not changing, monitored by observers whose detection slips by 1.5 per cent a year, produces an estimated trend of -0.0147 and a statistically significant decline in 39 per cent of programmes. A population that really is declining at three per cent, monitored the same way, gives -0.0449: the two trends simply add.

The estimator is doing its job. It is unbiased for the trend in the product, and the product is what was measured. No amount of design fixes this, because the design has no information about detection: the same annual means are consistent with a stable population and drifting observers, or a declining population and stable observers.

What does fix it is extra information, not extra effort. Repeat visits within a season support a detection model, a subset of sites can be surveyed by a second method for calibration, and observer identity recorded in the data at least allows the drift to be estimated rather than assumed away. All of these cost visits, which is why they belong in the design conversation rather than the analysis.

Where this leaves the series

The honest summary of these four posts is that a monitoring programme has two failure modes and they are asymmetric. The first is not seeing a real change, and that one is calculable in advance from the variance components: it is the power number, and it is mostly a question of how many years the programme survives. The second is reporting a change that is not there, and it does not come from the noise at all. It comes from the network changing under the programme, or the observers changing, and neither shows up in a standard error.

A design that names its detectable change, plans for attrition it can characterise, and records enough to estimate its own detection drift is a design that can be checked. One that promises to detect declines is not.

References

Legg and Nagy 2006 Journal of Environmental Management 78(2):194-199 (10.1016/j.jenvman.2005.04.016)

Nichols and Williams 2006 Trends in Ecology and Evolution 21(12):668-673 (10.1016/j.tree.2006.08.007)

Kery and Schmidt 2008 Community Ecology 9(2):207-216 (10.1556/ComEc.9.2008.2.10)

Gerrodette 1987 Ecology 68(5):1364-1372 (10.2307/1939220)

Related tutorials

• Variance components in monitoring data
• Revisit designs for long-term monitoring
• Power to detect a population trend
• Imperfect detection and occupancy