Checking a distance sampling model

distance sampling

model checking

abundance

ecology tutorial

The goodness-of-fit test accepts surveys whose density is 12 per cent wrong, and the assumption that matters most leaves no trace in the likelihood at all.

Author

Tidy Ecology

Published

2026-06-13

Three posts in this series have shown a distance sampling analysis returning a confident wrong number: a pooled detection function reversing which habitat is denser, a half-normal at a point inflating density by a sixth, a hierarchical model halving a population when the animals hold still.

The reasonable response is to ask what you would check. Distance sampling has a natural answer, because unlike most ecological models it makes a direct, testable prediction about data you can see: the histogram of distances should look like the fitted curve. So plot it, run a goodness-of-fit test, and if it passes, you are fine.

This post takes that programme seriously and measures how much protection it buys. The honest answer is: some, and a good deal less than the ritual implies.

library(ggplot2)
te_canvas <- theme(plot.background  = element_rect(fill = "#f5f4ee", colour = NA),
                   panel.background = element_rect(fill = "#f5f4ee", colour = NA))

g_hn  <- function(x, s) exp(-x^2/(2*s^2))
mu_hn <- function(s, w) sqrt(2*pi)*s*(pnorm(w, 0, s) - 0.5)
ld_hn <- function(x, s, w) -x^2/(2*s^2) - log(mu_hn(s, w))
g_hr  <- function(x, s, b) 1 - exp(-(x/s)^(-b))
mu_hr <- function(s, b, w) integrate(g_hr, 0, w, s = s, b = b)$value
fit_hn <- function(x, w) exp(optimise(function(p) -sum(ld_hn(x, exp(p), w)),
                                      c(-8, 1), tol = 1e-10)$minimum)
w <- 0.1

The check that works

Start with the good news, because there is some. Bin the distances, fit the key function, compare observed counts against expected counts with a chi-square statistic. Before trusting a test to reject a bad model, check that it accepts a good one at the advertised rate: a test with the wrong size is not evidence, it is a coin with an opinion.

cut <- seq(0, w, length.out = 6)
pi_hn <- function(s) { P <- sqrt(2*pi)*s*(pnorm(cut[-1], 0, s) - pnorm(cut[-6], 0, s))/w
                       P/sum(P) }
gof_hn <- function(x) {
  s  <- fit_hn(x, w)
  ob <- tabulate(findInterval(x, cut, rightmost.closed = TRUE), 5)
  ex <- length(x)*pi_hn(s)
  c(p = pchisq(sum((ob - ex)^2/ex), df = 5 - 1 - 1, lower.tail = FALSE), sigma = s)
}
## a world where the half-normal is exactly right
s_ok <- 0.045; D0 <- 100; L_ok <- 266/(2*D0*mu_hn(s_ok, w))
cal <- t(sapply(1:400, function(i) {
  set.seed(61000 + i)
  N <- rpois(1, D0*2*w*L_ok); x <- runif(N, 0, w)
  gof_hn(x[runif(N) < g_hn(x, s_ok)])
}))
sprintf("model-true data, n approx 266: test rejects at 0.05 in %.1f%% of surveys | median p %.3f",
        100*mean(cal[, "p"] < 0.05), median(cal[, "p"]))

[1] "model-true data, n approx 266: test rejects at 0.05 in 6.0% of surveys | median p 0.493"

Well behaved: it rejects the truth at about the nominal rate, so it is a real test rather than a formality. And it does catch gross problems. A spike of animals recorded at zero distance, a histogram with a hole in it, distances heaped on round numbers: all of that shows up.

What it misses

Now give it the job it is actually asked to do. The true detection function has a shoulder, generated by a hazard-rate key. You fit a half-normal, which cannot make a shoulder. From the point transect post, we know this costs about 11 per cent on a line. Does the test tell you?

s_hr <- 0.05; b_hr <- 3
mu_t <- mu_hr(s_hr, b_hr, w); Ln <- 266/(2*D0*mu_t)
res <- t(sapply(1:400, function(i) {
  set.seed(63000 + i)
  N <- rpois(1, D0*2*w*Ln); x <- runif(N, 0, w)
  x <- x[runif(N) < g_hr(x, s_hr, b_hr)]
  g <- gof_hn(x)
  c(p = g["p"], bias = length(x)/(2*Ln*mu_hn(g["sigma"], w))/D0 - 1, n = length(x))
}))
acc <- res[res[, 1] > 0.05, ]
sprintf("half-normal fitted to hazard-rate data (n approx %.0f, 5 bins, df 3):", mean(res[, 3]))

[1] "half-normal fitted to hazard-rate data (n approx 268, 5 bins, df 3):"

sprintf("  rejects at 0.05: %.1f%% of surveys | median p-value %.3f | mean density bias %+.1f%%",
        100*mean(res[, 1] < 0.05), median(res[, 1]), 100*mean(res[, 2]))

[1] "  rejects at 0.05: 17.2% of surveys | median p-value 0.312 | mean density bias +12.4%"

sprintf("  among the surveys the test ACCEPTS: mean bias %+.1f%%, worst %+.1f%%",
        100*mean(acc[, 2]), 100*max(acc[, 2]))

[1] "  among the surveys the test ACCEPTS: mean bias +12.2%, worst +37.8%"

sprintf("  correlation between the p-value and |bias|: %.3f", cor(res[, 1], abs(res[, 2])))

[1] "  correlation between the p-value and |bias|: 0.024"

Read those three lines slowly.

The test rejects the wrong model in 17 per cent of surveys, so 83 per cent of the time it hands the wrong model a clean bill of health. That is a power problem, and power problems are familiar and forgivable.

The second line is the one that should worry you. Among the surveys that pass, the mean density bias is +12.2 per cent, which is the same as the mean bias across all surveys. Passing the test tells you nothing about whether your density is right. The worst passing survey is out by 38 per cent, and it passed.

And the third line says it exactly: the correlation between the p-value and the size of the error is 0.024. Not weak. Zero, to within the noise of 400 surveys. The statistic the test computes and the quantity you care about are unrelated.

A goodness-of-fit test asks whether the fitted curve describes the shape of the histogram. Your density estimate depends on the height of that curve at zero. Those are different questions, and a model can answer the first one well while getting the second one wrong, because five bins of a smooth curve leave plenty of room to be wrong in the corner where nobody is looking.

The assumption with no fingerprints

That was the diagnostic being weak. This one is different in kind.

Every density estimate in this series divides by an effective width computed on the assumption that $g(0) = 1$: everything on the transect line is detected. It is the first assumption in every textbook, and it is the one that fails most often in practice. Divers miss animals below them. Aerial observers miss animals under the plane. A bird singing directly overhead in the canopy is not obviously easier to detect than one at 20 m.

Suppose it is false, and the true detection function is $c \cdot g(x)$ for some $c < 1$. Fit your model to the distances and see what happens.

set.seed(4283)
s_t <- 0.05; L <- 30; Dtrue <- 60
N <- rpois(1, Dtrue*2*w*L); xa <- runif(N, 0, w)
x <- xa[runif(N) < g_hn(xa, s_t)]
base <- optimise(function(p) -sum(ld_hn(x, exp(p), w)), c(-8, 1), tol = 1e-12)
out <- t(sapply(c(1, 0.8, 0.6, 0.4), function(g0) {
  ## the world where only a fraction g0 of the animals on the line are seen
  nll <- function(p) { s <- exp(p)
    -sum(log(g0*g_hn(x, s)) - log(g0*mu_hn(s, w))) }
  o <- optimise(nll, c(-8, 1), tol = 1e-12)
  c(g0 = g0, sigma = exp(o$minimum), d_sigma = exp(o$minimum) - exp(base$minimum),
    loglik = -o$objective, d_loglik = -o$objective + base$objective,
    D = length(x)/(2*L*mu_hn(exp(o$minimum), w)*g0))
}))
sprintf("one survey, n = %d detections", length(x))

[1] "one survey, n = 223 detections"

print(round(out, 8))

      g0      sigma d_sigma   loglik d_loglik         D
[1,] 1.0 0.05705162       0 531.7441        0  56.47631
[2,] 0.8 0.05705162       0 531.7441        0  70.59539
[3,] 0.6 0.05705162       0 531.7441        0  94.12719
[4,] 0.4 0.05705162       0 531.7441        0 141.19078

Look at the middle columns. The fitted scale is the same in all four rows, and so is the log-likelihood: those differences are not small, they are zero to the last bit an optimiser can carry. Meanwhile the density runs from 56.5 to 141.2, a factor of 2.5.

This is not a weak diagnostic. It is an exact invariance, and it takes one line of algebra to see why. The likelihood is built from the probability density of the observed distances,

\[f(x) = \frac{c\,g(x)}{\int_0^w c\,g(t)\,\mathrm{d}t} = \frac{g(x)}{\mu}\]

and the constant cancels. Your distances are conditional on detection, so they know the shape of the detection function and nothing whatsoever about its height. Every estimate in the series scales as $\hat{D} = D/c$: +25 per cent at $c = 0.8$, +67 per cent at $c = 0.6$, +150 per cent at $c = 0.4$.

No goodness-of-fit test can see this, because the fit is identical. No residual can see it. No AIC comparison can see it, because both models have the same likelihood and the same parameter count. There is no diagnostic, no plot and no statistic in R that will tell you, because the information is not in the data. The assumption you most need is the one your survey cannot check.

Left panel shows four detection curves at different heights, all producing one identical density of observed distances. Right panel shows the fitted sigma and log-likelihood flat across g(0) while density rises steeply. — Figure 2: Left: four detection functions differing only in their height at zero, and the single distribution of observed distances that all four produce. Right: what changes and what does not as g(0) falls.

What to check, and what to decide instead

Sorting the assumptions by whether the data can speak to them produces a short and lopsided table.

Assumption	Can the distances check it	What settles it
Key function has the right shape	Weakly	Goodness of fit, AIC, shouldered key by default
Distances measured accurately	Partly	Repeat measurement of a subset
No heaping on round numbers	Yes	The histogram, and a rangefinder
Detection depends on the recorded covariates	Partly	AIC, once the covariate is recorded
Animals uniform with respect to the line	No	Random placement of the design
g(0) = 1	No	Double observer or mark-recapture distance sampling
Animals do not move before detection	No	Snapshot protocol, or a faster observer

Two of the seven rows are about analysis. The other five are decisions made before anyone leaves the office, and the pattern is the one this blog keeps running into: the strongest thing you can do for an estimate is to make an assumption true by construction rather than to test it afterwards.

For $g(0)$ specifically, the fix is a second look. Two observers working independently, or a double platform, gives you a mark-recapture sample within the distance sampling survey: animals seen by one observer and missed by the other carry the information about what both of them missed. That information is not in a single observer’s distances at any sample size, and it never will be. Buckland et al. (2015) work through the designs; the point here is only that it is a field decision, not an analysis one.

The design is a field decision. What it buys is not, and that part is worth measuring rather than assuming. A second observer hands you a mark-recapture sample, and the sample needs an independence assumption of its own. Two are in use. Full independence treats the observers as independent at every distance. Point independence, which is the standard recommendation, assumes it only on the line and takes the shape of the curve from the pooled distances. Both are run below against the same surveys, in two worlds that carry the same amount of individual variation in detectability and differ only in whether that variation is still present at zero distance.

set.seed(90210)
do_g0 <- 0.8; do_s <- 0.05; do_N <- 1200; do_reps <- 150
do_grid <- seq(0, w, length.out = 400)

## match the two worlds on the spread of each animal's own detection probability
do_cv <- function(cc, ss) { p <- cc*mu_hn(ss, w)/w; sd(p)/mean(p) }
do_u   <- rnorm(2e5)
do_a0  <- uniroot(function(a) mean(plogis(a + 2*do_u)) - do_g0, c(-20, 20))$root
do_cvl <- do_cv(plogis(do_a0 + 2*do_u), do_s)
do_ts  <- uniroot(function(t) do_cv(do_g0, do_s*exp(t*do_u - t^2/2)) - do_cvl, c(0.01, 3))$root

do_sim <- function(at_line) {
  xx <- runif(do_N, 0, w)
  cc <- if (at_line) plogis(do_a0 + 2*rnorm(do_N)) else rep(do_g0, do_N)
  ss <- if (at_line) do_s else do_s*exp(do_ts*rnorm(do_N) - do_ts^2/2)
  p  <- cc*g_hn(xx, ss)              # one animal, one detectability, offered to both
  a  <- runif(do_N) < p; b <- runif(do_N) < p; k <- a | b
  list(x = xx[k], a = a[k], b = b[k])
}
do_key  <- function(xx, q) g_hn(xx, exp(q[1]))*(1 + q[2]*cos(pi*xx/w))/(1 + q[2])
do_fitk <- function(xx) {
  nk <- function(q) { gg <- do_key(do_grid, q); if (any(gg <= 0)) return(1e9)
    -(sum(log(pmax(do_key(xx, q), 1e-12))) - length(xx)*log(mean(gg)*w)) }
  o <- optim(c(log(do_s), 0.1), nk); optim(o$par, nk)$par
}
do_cds <- function(xx) length(xx)*w/mu_hn(fit_hn(xx, w), w)     # one observer, g(0) = 1
do_full <- function(d) {                                        # MRDS, full independence
  nll <- function(q) { th <- plogis(q[1]); s <- exp(q[2]); p <- th*g_hn(d$x, s)
    Pd <- (2*th*mu_hn(s, w) - th^2*mu_hn(s/sqrt(2), w))/w
    -(sum(ifelse(d$a & d$b, 2*log(p), log(p) + log1p(-p))) - length(d$x)*log(Pd)) }
  o <- optim(c(1, log(do_s)), nll); o <- optim(o$par, nll)
  sum(1/(1 - (1 - plogis(o$par[1])*g_hn(d$x, exp(o$par[2])))^2))
}
do_p0 <- function(d) {                    # height, from the mark-recapture model at zero
  xr <- c(d$x[d$a], d$x[d$b]); yr <- c(d$b[d$a], d$a[d$b])
  m <- optim(c(1, -1, 0), function(bb) { e <- bb[1] + bb[2]*xr/w + bb[3]*(xr/w)^2
                                         -sum(yr*e - log1p(exp(e))) })
  1 - (1 - plogis(m$par[1]))^2
}
do_point <- function(d, adjust) {                 # MRDS, point independence
  pd0 <- do_p0(d)
  if (adjust) { q <- do_fitk(d$x); gg <- do_key(d$x, q)/do_key(0, q) }
  else gg <- g_hn(d$x, fit_hn(d$x, w))
  sum(1/(pd0*gg))
}
do_run <- function(at_line) t(replicate(do_reps, { d <- do_sim(at_line)
  c(conv = do_cds(d$x[d$a])/do_N, full = do_full(d)/do_N,
    point = do_point(d, TRUE)/do_N, bare = do_point(d, FALSE)/do_N) }))
do_L <- colMeans(do_run(TRUE)); do_S <- colMeans(do_run(FALSE))
do_rec <- function(v, k) unname((v[k] - v["conv"])/(1 - v["conv"]))
do_rec_line_full   <- do_rec(do_L, "full");  do_rec_line_point  <- do_rec(do_L, "point")
do_rec_shape_full  <- do_rec(do_S, "full");  do_rec_shape_point <- do_rec(do_S, "point")
do_rec_line_bare   <- do_rec(do_L, "bare");  do_rec_shape_bare  <- do_rec(do_S, "bare")

## pooling robustness on its own terms: g(0) = 1 genuinely true, same spread of scale
do_pool <- function(adjust) mean(replicate(do_reps, {
  xx <- runif(do_N, 0, w); ss <- do_s*exp(do_ts*rnorm(do_N) - do_ts^2/2)
  dd <- xx[runif(do_N) < g_hn(xx, ss)]
  if (adjust) { q <- do_fitk(dd); gg <- do_key(do_grid, q)/do_key(0, q)
                length(dd)/(mean(gg)*w)*w/do_N } else do_cds(dd)/do_N }))
do_pool_bare <- do_pool(FALSE); do_pool_adj <- do_pool(TRUE)

cat(sprintf("both worlds carry the same CV of individual detection probability: %.2f\n", do_cvl))

both worlds carry the same CV of individual detection probability: 0.29

for (i in 1:2) { v <- list(do_L, do_S)[[i]]
  cat(sprintf("\n%s\n", c("heterogeneity SURVIVING at zero distance",
                          "heterogeneity that VANISHES at zero distance")[i]))
  cat(sprintf("  one observer, g(0) = 1 assumed    : N_hat/N %.3f\n", v["conv"]))
  cat(sprintf("  two observers, full independence  : N_hat/N %.3f, gap recovered %.2f\n",
              v["full"], do_rec(v, "full")))
  cat(sprintf("  two observers, point independence : N_hat/N %.3f, gap recovered %.2f\n",
              v["point"], do_rec(v, "point")))
  cat(sprintf("     the same, pooled key left bare : N_hat/N %.3f, gap recovered %.2f\n",
              v["bare"], do_rec(v, "bare"))) }


heterogeneity SURVIVING at zero distance
  one observer, g(0) = 1 assumed    : N_hat/N 0.803
  two observers, full independence  : N_hat/N 0.919, gap recovered 0.59
  two observers, point independence : N_hat/N 0.913, gap recovered 0.56
     the same, pooled key left bare : N_hat/N 0.966, gap recovered 0.83

heterogeneity that VANISHES at zero distance
  one observer, g(0) = 1 assumed    : N_hat/N 0.747
  two observers, full independence  : N_hat/N 0.891, gap recovered 0.57
  two observers, point independence : N_hat/N 0.980, gap recovered 0.92
     the same, pooled key left bare : N_hat/N 0.981, gap recovered 0.92

cat(sprintf("\npooling robustness alone, g(0) = 1 true: bare key %.3f, adjusted key %.3f\n",
            do_pool_bare, do_pool_adj))


pooling robustness alone, g(0) = 1 true: bare key 0.935, adjusted key 0.952

Same mean $g(0)$, same amount of heterogeneity, two different answers. Where the variation survives on the line, the two estimators agree and both fall short: point independence recovers 0.56 of the gap and full independence 0.59. Where it lives only in how far each animal can be seen, so that every animal on the line is equally detectable, point independence recovers 0.92 while full independence stays at 0.57. What decides the answer is where the heterogeneity lives, not how much of it there is, and that is a fact about your species and your observers rather than one your distances will settle.

Two of those statements are theorems and should not be sold as findings. That unmodelled individual heterogeneity pulls abundance down is the model-Mh direction, which the hierarchical distance sampling post sets out; it arrives here at zero distance rather than between visits. That point independence is exactly right when the heterogeneity vanishes on the line follows from its own definition, because the assumption it makes is then true, and full independence is exactly right in the mirror case where each observer’s heterogeneity is private to that observer. What is measured above is neither of those: it is how much of the gap each one closes when the conditions hold only partly, and how far apart the two get.

One warning goes with the point independence estimator. The shape of the curve it uses comes from the pooled distances, and two observers together always give a shouldered curve even when neither of them does alone. Fit that with a bare half-normal and the recovered fraction in the first world rises from 0.56 to 0.83, which is not the estimator getting better. A key that cannot make a shoulder puts the missing shoulder into the height, the upward error this post measured in its own second section at +12.4 per cent, and here it happens to point against the dependence bias. Report the second number and you have reported the cancellation of two mistakes.

The honest close

Distance sampling earns its reputation. It needs one visit rather than three. It gives an absolute density rather than an index. Its central property is real, with the conditions its own source attaches: pooling unmodelled heterogeneity leaves the population total close to unbiased, extreme heterogeneity excepted, and it does not extend to subpopulation estimates read off a single pooled detection function (Rexstad et al. 2023). With $g(0) = 1$ genuinely true and the spread of detectability used above, a bare half-normal returns 93.47 per cent of the population, and a cosine adjustment free to follow the pooled shape returns 95.24 per cent. The flexible key closes 1.77 of those 6.53 percentage points, so at a detection CV of 0.29 what stands between this estimate and unbiasedness is the amount of heterogeneity rather than the shape of the key. Close to unbiased is a genuinely unusual property for an estimator to have, and it is not the same as unbiased.

The failures in this series are consistent with all of that, and they share a shape. The estimator reads one number off a curve at the origin, and everything that goes wrong goes wrong by moving that number: pooling moves it the same way for every stratum and destroys the comparison while sparing the total; the wrong key moves it more at a point than on a line; a missing $g(0)$ scales it by a constant nobody can estimate from distances. In every case the fit stays clean, the standard error stays narrow, and the output looks exactly like an output that is right.

Which is the sentence this cluster has been building towards, and it is not really about distance sampling: the diagnostics available to you are the ones that operate on the data you collected, and the assumptions that hurt you most are the ones that determined what data you were able to collect. Check what you can. Decide the rest with a protocol, write down which is which, and report the number with the assumption attached rather than with an interval that has quietly forgotten it.

References

Buckland ST 2006 The Auk 123(2):345-357 (10.1093/auk/123.2.345)

Buckland ST, Rexstad EA, Marques TA, Oedekoven CS 2015 Distance Sampling: Methods and Applications. Springer. ISBN 978-3-319-19218-5

Marques FFC, Buckland ST 2003 Biometrics 59(4):924-935 (10.1111/j.0006-341X.2003.00107.x)

Rexstad E, Buckland ST, Marshall L, Borchers D 2023 Ecology and Evolution 13(1):e9684 (10.1002/ece3.9684)

--- title: "Checking a distance sampling model" description: "The goodness-of-fit test accepts surveys whose density is 12 per cent wrong, and the assumption that matters most leaves no trace in the likelihood at all." date: "2026-06-13 10:00" categories: [R, distance sampling, model checking, abundance, ecology tutorial] image: thumbnail.png image-alt: "A scatter of goodness-of-fit p-values against density bias showing no relationship" --- Three posts in this series have shown a distance sampling analysis returning a confident wrong number: [a pooled detection function](../covariates-in-the-detection-function/) reversing which habitat is denser, [a half-normal at a point](../point-transect-distance-sampling/) inflating density by a sixth, [a hierarchical model](../hierarchical-distance-sampling/) halving a population when the animals hold still. The reasonable response is to ask what you would check. Distance sampling has a natural answer, because unlike most ecological models it makes a direct, testable prediction about data you can see: the histogram of distances should look like the fitted curve. So plot it, run a goodness-of-fit test, and if it passes, you are fine. This post takes that programme seriously and measures how much protection it buys. The honest answer is: some, and a good deal less than the ritual implies. ```{r setup} #| message: false #| warning: false library(ggplot2) te_canvas <- theme(plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA)) g_hn <- function(x, s) exp(-x^2/(2*s^2)) mu_hn <- function(s, w) sqrt(2*pi)*s*(pnorm(w, 0, s) - 0.5) ld_hn <- function(x, s, w) -x^2/(2*s^2) - log(mu_hn(s, w)) g_hr <- function(x, s, b) 1 - exp(-(x/s)^(-b)) mu_hr <- function(s, b, w) integrate(g_hr, 0, w, s = s, b = b)$value fit_hn <- function(x, w) exp(optimise(function(p) -sum(ld_hn(x, exp(p), w)), c(-8, 1), tol = 1e-10)$minimum) w <- 0.1 ``` ## The check that works Start with the good news, because there is some. Bin the distances, fit the key function, compare observed counts against expected counts with a chi-square statistic. Before trusting a test to reject a bad model, check that it accepts a good one at the advertised rate: a test with the wrong size is not evidence, it is a coin with an opinion. ```{r calibrate} cut <- seq(0, w, length.out = 6) pi_hn <- function(s) { P <- sqrt(2*pi)*s*(pnorm(cut[-1], 0, s) - pnorm(cut[-6], 0, s))/w P/sum(P) } gof_hn <- function(x) { s <- fit_hn(x, w) ob <- tabulate(findInterval(x, cut, rightmost.closed = TRUE), 5) ex <- length(x)*pi_hn(s) c(p = pchisq(sum((ob - ex)^2/ex), df = 5 - 1 - 1, lower.tail = FALSE), sigma = s) } ## a world where the half-normal is exactly right s_ok <- 0.045; D0 <- 100; L_ok <- 266/(2*D0*mu_hn(s_ok, w)) cal <- t(sapply(1:400, function(i) { set.seed(61000 + i) N <- rpois(1, D0*2*w*L_ok); x <- runif(N, 0, w) gof_hn(x[runif(N) < g_hn(x, s_ok)]) })) sprintf("model-true data, n approx 266: test rejects at 0.05 in %.1f%% of surveys | median p %.3f", 100*mean(cal[, "p"] < 0.05), median(cal[, "p"])) ``` Well behaved: it rejects the truth at about the nominal rate, so it is a real test rather than a formality. And it does catch gross problems. A spike of animals recorded at zero distance, a histogram with a hole in it, distances heaped on round numbers: all of that shows up. ## What it misses Now give it the job it is actually asked to do. The true detection function has a shoulder, generated by a hazard-rate key. You fit a half-normal, which cannot make a shoulder. From the [point transect post](../point-transect-distance-sampling/), we know this costs about 11 per cent on a line. Does the test tell you? ```{r blind} s_hr <- 0.05; b_hr <- 3 mu_t <- mu_hr(s_hr, b_hr, w); Ln <- 266/(2*D0*mu_t) res <- t(sapply(1:400, function(i) { set.seed(63000 + i) N <- rpois(1, D0*2*w*Ln); x <- runif(N, 0, w) x <- x[runif(N) < g_hr(x, s_hr, b_hr)] g <- gof_hn(x) c(p = g["p"], bias = length(x)/(2*Ln*mu_hn(g["sigma"], w))/D0 - 1, n = length(x)) })) acc <- res[res[, 1] > 0.05, ] sprintf("half-normal fitted to hazard-rate data (n approx %.0f, 5 bins, df 3):", mean(res[, 3])) sprintf(" rejects at 0.05: %.1f%% of surveys | median p-value %.3f | mean density bias %+.1f%%", 100*mean(res[, 1] < 0.05), median(res[, 1]), 100*mean(res[, 2])) sprintf(" among the surveys the test ACCEPTS: mean bias %+.1f%%, worst %+.1f%%", 100*mean(acc[, 2]), 100*max(acc[, 2])) sprintf(" correlation between the p-value and |bias|: %.3f", cor(res[, 1], abs(res[, 2]))) ``` Read those three lines slowly. The test rejects the wrong model in `r sprintf("%.0f", 100*mean(res[, 1] < 0.05))` per cent of surveys, so `r sprintf("%.0f", 100*mean(res[, 1] > 0.05))` per cent of the time it hands the wrong model a clean bill of health. That is a power problem, and power problems are familiar and forgivable. The second line is the one that should worry you. Among the surveys that *pass*, the mean density bias is `r sprintf("%+.1f", 100*mean(acc[, 2]))` per cent, which is the same as the mean bias across all surveys. Passing the test tells you nothing about whether your density is right. The worst passing survey is out by `r sprintf("%.0f", 100*max(acc[, 2]))` per cent, and it passed. And the third line says it exactly: the correlation between the p-value and the size of the error is `r sprintf("%.3f", cor(res[, 1], abs(res[, 2])))`. Not weak. Zero, to within the noise of 400 surveys. The statistic the test computes and the quantity you care about are unrelated. > A goodness-of-fit test asks whether the fitted curve describes the shape of the histogram. Your density estimate depends on the height of that curve at zero. Those are different questions, and a model can answer the first one well while getting the second one wrong, because five bins of a smooth curve leave plenty of room to be wrong in the corner where nobody is looking. ```{r fig-blind} #| echo: false #| fig-width: 8 #| fig-height: 3.6 #| fig-cap: "Four hundred surveys of a shouldered detection function analysed with a half-normal. Left: the p-value against the resulting density error, with the accepted surveys in green. Right: the distribution of density error among surveys the test accepts, against the truth." #| fig-alt: "Left panel is a scatter of goodness-of-fit p-values against density bias showing a formless cloud with no trend. Right panel shows a histogram of density bias among accepted surveys, bulk of it well above the true value with a thin tail reaching below." ink <- "#16241d"; forest <- "#275139"; gold <- "#c9b458"; rust <- "#b5534e" dd <- data.frame(p = res[, 1], bias = 100*res[, 2], k = ifelse(res[, 1] > 0.05, "test accepts", "test rejects")) p1 <- ggplot(dd, aes(p, bias, colour = k)) + geom_hline(yintercept = 0, colour = ink, linetype = 2) + geom_vline(xintercept = 0.05, colour = rust, linetype = 3) + geom_point(size = 1, alpha = 0.7) + scale_colour_manual(values = c("test accepts" = forest, "test rejects" = rust)) + labs(x = "goodness-of-fit p-value", y = "density error (%)", colour = NULL, title = "No relationship at all", subtitle = sprintf("correlation of p with |error|: %.3f", cor(res[, 1], abs(res[, 2])))) + theme_minimal(base_size = 10) + te_canvas + theme(legend.position = "bottom", plot.title = element_text(colour = ink, face = "bold"), plot.subtitle = element_text(size = 8)) p2 <- ggplot(dd[dd$k == "test accepts", ], aes(bias)) + geom_histogram(bins = 30, fill = forest, alpha = 0.6) + geom_vline(xintercept = 0, colour = ink, linetype = 2, linewidth = 0.7) + labs(x = "density error (%)", y = "surveys", title = "The surveys that passed", subtitle = sprintf("mean error %+.1f%%, only %.0f%% within 2%% of the truth", 100*mean(acc[, 2]), 100*mean(abs(acc[, 2]) < 0.02))) + theme_minimal(base_size = 10) + te_canvas + theme(plot.title = element_text(colour = ink, face = "bold"), plot.subtitle = element_text(size = 8)) 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)) ``` ## The assumption with no fingerprints That was the diagnostic being weak. This one is different in kind. Every density estimate in this series divides by an effective width computed on the assumption that $g(0) = 1$: everything on the transect line is detected. It is the first assumption in every textbook, and it is the one that fails most often in practice. Divers miss animals below them. Aerial observers miss animals under the plane. A bird singing directly overhead in the canopy is not obviously easier to detect than one at 20 m. Suppose it is false, and the true detection function is $c \cdot g(x)$ for some $c < 1$. Fit your model to the distances and see what happens. ```{r g0} set.seed(4283) s_t <- 0.05; L <- 30; Dtrue <- 60 N <- rpois(1, Dtrue*2*w*L); xa <- runif(N, 0, w) x <- xa[runif(N) < g_hn(xa, s_t)] base <- optimise(function(p) -sum(ld_hn(x, exp(p), w)), c(-8, 1), tol = 1e-12) out <- t(sapply(c(1, 0.8, 0.6, 0.4), function(g0) { ## the world where only a fraction g0 of the animals on the line are seen nll <- function(p) { s <- exp(p) -sum(log(g0*g_hn(x, s)) - log(g0*mu_hn(s, w))) } o <- optimise(nll, c(-8, 1), tol = 1e-12) c(g0 = g0, sigma = exp(o$minimum), d_sigma = exp(o$minimum) - exp(base$minimum), loglik = -o$objective, d_loglik = -o$objective + base$objective, D = length(x)/(2*L*mu_hn(exp(o$minimum), w)*g0)) })) sprintf("one survey, n = %d detections", length(x)) print(round(out, 8)) ``` Look at the middle columns. The fitted scale is the same in all four rows, and so is the log-likelihood: those differences are not small, they are zero to the last bit an optimiser can carry. Meanwhile the density runs from `r sprintf("%.1f", out[1, "D"])` to `r sprintf("%.1f", out[4, "D"])`, a factor of `r sprintf("%.1f", out[4, "D"]/out[1, "D"])`. This is not a weak diagnostic. It is an exact invariance, and it takes one line of algebra to see why. The likelihood is built from the probability density of the *observed* distances, $$f(x) = \frac{c\,g(x)}{\int_0^w c\,g(t)\,\mathrm{d}t} = \frac{g(x)}{\mu}$$ and the constant cancels. Your distances are conditional on detection, so they know the *shape* of the detection function and nothing whatsoever about its height. Every estimate in the series scales as $\hat{D} = D/c$: `r sprintf("%+.0f", 100*(1/0.8 - 1))` per cent at $c = 0.8$, `r sprintf("%+.0f", 100*(1/0.6 - 1))` per cent at $c = 0.6$, `r sprintf("%+.0f", 100*(1/0.4 - 1))` per cent at $c = 0.4$. No goodness-of-fit test can see this, because the fit is identical. No residual can see it. No AIC comparison can see it, because both models have the same likelihood and the same parameter count. There is no diagnostic, no plot and no statistic in R that will tell you, because the information is not in the data. The assumption you most need is the one your survey cannot check. ```{r fig-g0} #| echo: false #| fig-width: 8 #| fig-height: 3.6 #| fig-cap: "Left: four detection functions differing only in their height at zero, and the single distribution of observed distances that all four produce. Right: what changes and what does not as g(0) falls." #| fig-alt: "Left panel shows four detection curves at different heights, all producing one identical density of observed distances. Right panel shows the fitted sigma and log-likelihood flat across g(0) while density rises steeply." xs <- seq(0, w, length.out = 300) cur <- do.call(rbind, lapply(c(1, 0.8, 0.6, 0.4), function(g0) data.frame(x = xs, g = g0*g_hn(xs, s_t), k = sprintf("g(0) = %.1f", g0)))) p1 <- ggplot(cur, aes(x*1000, g, colour = k)) + geom_line(linewidth = 0.9) + geom_line(data = data.frame(x = xs, g = g_hn(xs, s_t)/mu_hn(s_t, w)/20.2, k = "observed f(x)"), colour = ink, linetype = 2, linewidth = 1) + scale_colour_manual(values = c("g(0) = 1.0" = forest, "g(0) = 0.8" = gold, "g(0) = 0.6" = rust, "g(0) = 0.4" = "#93a87f")) + labs(x = "perpendicular distance (m)", y = "detection probability", colour = NULL, title = "Four worlds, one dataset", subtitle = "dashed: the distances you record, scaled; identical in all four") + theme_minimal(base_size = 10) + te_canvas + theme(legend.position = "bottom", plot.title = element_text(colour = ink, face = "bold"), plot.subtitle = element_text(size = 8), legend.text = element_text(size = 7)) sc <- data.frame(g0 = rep(out[, "g0"], 3), v = c(out[, "sigma"]/out[1, "sigma"], out[, "loglik"]/out[1, "loglik"], out[, "D"]/out[1, "D"]), k = rep(c("fitted sigma", "log-likelihood", "estimated density"), each = 4)) p2 <- ggplot(sc, aes(g0, v, colour = k)) + geom_hline(yintercept = 1, colour = ink, linetype = 3) + geom_line(linewidth = 0.9) + geom_point(size = 2) + scale_colour_manual(values = c("fitted sigma" = forest, "log-likelihood" = gold, "estimated density" = rust)) + scale_x_reverse() + labs(x = "true g(0)", y = "relative to the g(0) = 1 analysis", colour = NULL, title = "Everything you can check is flat", subtitle = "the one thing you report is not") + theme_minimal(base_size = 10) + te_canvas + theme(legend.position = "bottom", plot.title = element_text(colour = ink, face = "bold"), plot.subtitle = element_text(size = 8)) 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 to check, and what to decide instead Sorting the assumptions by whether the data can speak to them produces a short and lopsided table. ```{r checklist} #| echo: false knitr::kable(data.frame( Assumption = c("Key function has the right shape", "Distances measured accurately", "No heaping on round numbers", "Detection depends on the recorded covariates", "Animals uniform with respect to the line", "g(0) = 1", "Animals do not move before detection"), `Can the distances check it` = c("Weakly", "Partly", "Yes", "Partly", "No", "No", "No"), `What settles it` = c("Goodness of fit, AIC, shouldered key by default", "Repeat measurement of a subset", "The histogram, and a rangefinder", "AIC, once the covariate is recorded", "Random placement of the design", "Double observer or mark-recapture distance sampling", "Snapshot protocol, or a faster observer"), check.names = FALSE)) ``` Two of the seven rows are about analysis. The other five are decisions made before anyone leaves the office, and the pattern is the one this blog keeps running into: the strongest thing you can do for an estimate is to make an assumption true by construction rather than to test it afterwards. For $g(0)$ specifically, the fix is a second look. Two observers working independently, or a double platform, gives you a mark-recapture sample within the distance sampling survey: animals seen by one observer and missed by the other carry the information about what both of them missed. That information is not in a single observer's distances at any sample size, and it never will be. Buckland et al. (2015) work through the designs; the point here is only that it is a field decision, not an analysis one. The design is a field decision. What it buys is not, and that part is worth measuring rather than assuming. A second observer hands you a mark-recapture sample, and the sample needs an independence assumption of its own. Two are in use. Full independence treats the observers as independent at every distance. Point independence, which is the standard recommendation, assumes it only on the line and takes the shape of the curve from the pooled distances. Both are run below against the same surveys, in two worlds that carry the same amount of individual variation in detectability and differ only in whether that variation is still present at zero distance. ```{r double-observer} set.seed(90210) do_g0 <- 0.8; do_s <- 0.05; do_N <- 1200; do_reps <- 150 do_grid <- seq(0, w, length.out = 400) ## match the two worlds on the spread of each animal's own detection probability do_cv <- function(cc, ss) { p <- cc*mu_hn(ss, w)/w; sd(p)/mean(p) } do_u <- rnorm(2e5) do_a0 <- uniroot(function(a) mean(plogis(a + 2*do_u)) - do_g0, c(-20, 20))$root do_cvl <- do_cv(plogis(do_a0 + 2*do_u), do_s) do_ts <- uniroot(function(t) do_cv(do_g0, do_s*exp(t*do_u - t^2/2)) - do_cvl, c(0.01, 3))$root do_sim <- function(at_line) { xx <- runif(do_N, 0, w) cc <- if (at_line) plogis(do_a0 + 2*rnorm(do_N)) else rep(do_g0, do_N) ss <- if (at_line) do_s else do_s*exp(do_ts*rnorm(do_N) - do_ts^2/2) p <- cc*g_hn(xx, ss) # one animal, one detectability, offered to both a <- runif(do_N) < p; b <- runif(do_N) < p; k <- a | b list(x = xx[k], a = a[k], b = b[k]) } do_key <- function(xx, q) g_hn(xx, exp(q[1]))*(1 + q[2]*cos(pi*xx/w))/(1 + q[2]) do_fitk <- function(xx) { nk <- function(q) { gg <- do_key(do_grid, q); if (any(gg <= 0)) return(1e9) -(sum(log(pmax(do_key(xx, q), 1e-12))) - length(xx)*log(mean(gg)*w)) } o <- optim(c(log(do_s), 0.1), nk); optim(o$par, nk)$par } do_cds <- function(xx) length(xx)*w/mu_hn(fit_hn(xx, w), w) # one observer, g(0) = 1 do_full <- function(d) { # MRDS, full independence nll <- function(q) { th <- plogis(q[1]); s <- exp(q[2]); p <- th*g_hn(d$x, s) Pd <- (2*th*mu_hn(s, w) - th^2*mu_hn(s/sqrt(2), w))/w -(sum(ifelse(d$a & d$b, 2*log(p), log(p) + log1p(-p))) - length(d$x)*log(Pd)) } o <- optim(c(1, log(do_s)), nll); o <- optim(o$par, nll) sum(1/(1 - (1 - plogis(o$par[1])*g_hn(d$x, exp(o$par[2])))^2)) } do_p0 <- function(d) { # height, from the mark-recapture model at zero xr <- c(d$x[d$a], d$x[d$b]); yr <- c(d$b[d$a], d$a[d$b]) m <- optim(c(1, -1, 0), function(bb) { e <- bb[1] + bb[2]*xr/w + bb[3]*(xr/w)^2 -sum(yr*e - log1p(exp(e))) }) 1 - (1 - plogis(m$par[1]))^2 } do_point <- function(d, adjust) { # MRDS, point independence pd0 <- do_p0(d) if (adjust) { q <- do_fitk(d$x); gg <- do_key(d$x, q)/do_key(0, q) } else gg <- g_hn(d$x, fit_hn(d$x, w)) sum(1/(pd0*gg)) } do_run <- function(at_line) t(replicate(do_reps, { d <- do_sim(at_line) c(conv = do_cds(d$x[d$a])/do_N, full = do_full(d)/do_N, point = do_point(d, TRUE)/do_N, bare = do_point(d, FALSE)/do_N) })) do_L <- colMeans(do_run(TRUE)); do_S <- colMeans(do_run(FALSE)) do_rec <- function(v, k) unname((v[k] - v["conv"])/(1 - v["conv"])) do_rec_line_full <- do_rec(do_L, "full"); do_rec_line_point <- do_rec(do_L, "point") do_rec_shape_full <- do_rec(do_S, "full"); do_rec_shape_point <- do_rec(do_S, "point") do_rec_line_bare <- do_rec(do_L, "bare"); do_rec_shape_bare <- do_rec(do_S, "bare") ## pooling robustness on its own terms: g(0) = 1 genuinely true, same spread of scale do_pool <- function(adjust) mean(replicate(do_reps, { xx <- runif(do_N, 0, w); ss <- do_s*exp(do_ts*rnorm(do_N) - do_ts^2/2) dd <- xx[runif(do_N) < g_hn(xx, ss)] if (adjust) { q <- do_fitk(dd); gg <- do_key(do_grid, q)/do_key(0, q) length(dd)/(mean(gg)*w)*w/do_N } else do_cds(dd)/do_N })) do_pool_bare <- do_pool(FALSE); do_pool_adj <- do_pool(TRUE) cat(sprintf("both worlds carry the same CV of individual detection probability: %.2f\n", do_cvl)) for (i in 1:2) { v <- list(do_L, do_S)[[i]] cat(sprintf("\n%s\n", c("heterogeneity SURVIVING at zero distance", "heterogeneity that VANISHES at zero distance")[i])) cat(sprintf(" one observer, g(0) = 1 assumed : N_hat/N %.3f\n", v["conv"])) cat(sprintf(" two observers, full independence : N_hat/N %.3f, gap recovered %.2f\n", v["full"], do_rec(v, "full"))) cat(sprintf(" two observers, point independence : N_hat/N %.3f, gap recovered %.2f\n", v["point"], do_rec(v, "point"))) cat(sprintf(" the same, pooled key left bare : N_hat/N %.3f, gap recovered %.2f\n", v["bare"], do_rec(v, "bare"))) } cat(sprintf("\npooling robustness alone, g(0) = 1 true: bare key %.3f, adjusted key %.3f\n", do_pool_bare, do_pool_adj)) ``` Same mean $g(0)$, same amount of heterogeneity, two different answers. Where the variation survives on the line, the two estimators agree and both fall short: point independence recovers `r sprintf("%.2f", do_rec_line_point)` of the gap and full independence `r sprintf("%.2f", do_rec_line_full)`. Where it lives only in how far each animal can be seen, so that every animal on the line is equally detectable, point independence recovers `r sprintf("%.2f", do_rec_shape_point)` while full independence stays at `r sprintf("%.2f", do_rec_shape_full)`. What decides the answer is where the heterogeneity lives, not how much of it there is, and that is a fact about your species and your observers rather than one your distances will settle. Two of those statements are theorems and should not be sold as findings. That unmodelled individual heterogeneity pulls abundance down is the model-Mh direction, which the [hierarchical distance sampling post](../hierarchical-distance-sampling/) sets out; it arrives here at zero distance rather than between visits. That point independence is exactly right when the heterogeneity vanishes on the line follows from its own definition, because the assumption it makes is then true, and full independence is exactly right in the mirror case where each observer's heterogeneity is private to that observer. What is measured above is neither of those: it is how much of the gap each one closes when the conditions hold only partly, and how far apart the two get. One warning goes with the point independence estimator. The shape of the curve it uses comes from the pooled distances, and two observers together always give a shouldered curve even when neither of them does alone. Fit that with a bare half-normal and the recovered fraction in the first world rises from `r sprintf("%.2f", do_rec_line_point)` to `r sprintf("%.2f", do_rec_line_bare)`, which is not the estimator getting better. A key that cannot make a shoulder puts the missing shoulder into the height, the upward error this post measured in its own second section at `r sprintf("%+.1f", 100*mean(res[, 2]))` per cent, and here it happens to point against the dependence bias. Report the second number and you have reported the cancellation of two mistakes. ## The honest close Distance sampling earns its reputation. It needs one visit rather than three. It gives an absolute density rather than an index. Its central property is real, with the conditions its own source attaches: pooling unmodelled heterogeneity leaves the population total close to unbiased, extreme heterogeneity excepted, and it does not extend to subpopulation estimates read off a single pooled detection function (Rexstad et al. 2023). With $g(0) = 1$ genuinely true and the spread of detectability used above, a bare half-normal returns `r sprintf("%.2f", 100*do_pool_bare)` per cent of the population, and a cosine adjustment free to follow the pooled shape returns `r sprintf("%.2f", 100*do_pool_adj)` per cent. The flexible key closes `r sprintf("%.2f", 100*(do_pool_adj - do_pool_bare))` of those `r sprintf("%.2f", 100*(1 - do_pool_bare))` percentage points, so at a detection CV of `r sprintf("%.2f", do_cvl)` what stands between this estimate and unbiasedness is the amount of heterogeneity rather than the shape of the key. Close to unbiased is a genuinely unusual property for an estimator to have, and it is not the same as unbiased. The failures in this series are consistent with all of that, and they share a shape. The estimator reads one number off a curve at the origin, and everything that goes wrong goes wrong by moving that number: pooling moves it the same way for every stratum and destroys the comparison while sparing the total; the wrong key moves it more at a point than on a line; a missing $g(0)$ scales it by a constant nobody can estimate from distances. In every case the fit stays clean, the standard error stays narrow, and the output looks exactly like an output that is right. Which is the sentence this cluster has been building towards, and it is not really about distance sampling: the diagnostics available to you are the ones that operate on the data you collected, and the assumptions that hurt you most are the ones that determined what data you were able to collect. Check what you can. Decide the rest with a protocol, write down which is which, and report the number with the assumption attached rather than with an interval that has quietly forgotten it. ## References Buckland ST 2006 The Auk 123(2):345-357 (10.1093/auk/123.2.345) Buckland ST, Rexstad EA, Marques TA, Oedekoven CS 2015 Distance Sampling: Methods and Applications. Springer. ISBN 978-3-319-19218-5 Marques FFC, Buckland ST 2003 Biometrics 59(4):924-935 (10.1111/j.0006-341X.2003.00107.x) Rexstad E, Buckland ST, Marshall L, Borchers D 2023 Ecology and Evolution 13(1):e9684 (10.1002/ece3.9684) ## Related tutorials - [Distance sampling on line transects](../distance-sampling-line-transect/) - [Covariates in the detection function](../covariates-in-the-detection-function/) - [Digit preference and heaping in field data](../digit-preference-and-heaping/) - [Hierarchical distance sampling](../hierarchical-distance-sampling/)