---
title: "Hierarchical distance sampling"
description: "One visit with distances beats three visits of counts, and the reason is a biological bet about whether your animals hold still between visits."
date: 2026-08-26 11:00
categories: [R, distance sampling, abundance, imperfect detection, ecology tutorial]
image: thumbnail.png
image-alt: "Binned distance counts across survey sites with a fitted hierarchical model overlaid"
---
Two families of models estimate abundance from counts that miss animals, and they buy their detection estimate in different currencies. [N-mixture models](../n-mixture-models-abundance/) buy it with repeat visits: go back three times, and the variation between visits tells you how often you miss things. [Distance sampling](../distance-sampling-line-transect/) buys it with a ruler: the shape of the distance distribution tells you the same thing from a single pass.
Hierarchical distance sampling (Royle, Dawson and Bates 2004) is the second currency in the first family's shape. You get a site-level abundance model with covariates, exactly as in an N-mixture, and the detection parameter comes from the distances rather than from going back. This post builds it in base R, then runs the two designs against each other in a world where the assumptions differ, which turns out to be the whole story.
```{r setup}
#| message: false
#| warning: false
library(ggplot2)
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)
w <- 0.1 # truncation, km
R <- 200 # transect segments (sites)
ell <- 0.5 # segment length, km
A <- 2*w*ell # covered area per site, km2
cut <- seq(0, w, length.out = 6); B <- 5 # five 20 m distance bins
## cell probability: the chance a site's animal is detected AND lands in bin b
pi_b <- function(s) sqrt(2*pi)*s*(pnorm(cut[-1], 0, s) - pnorm(cut[-length(cut)], 0, s))/w
sprintf("sigma 0.05: bin probabilities %s | missed %.4f | detected %.4f (= mu/w = %.4f)",
paste(sprintf("%.4f", pi_b(0.05)), collapse = " "),
1 - sum(pi_b(0.05)), sum(pi_b(0.05)), mu_hn(0.05, w)/w)
```
## The identity that makes it easy
The model has a latent variable in it. Site $s$ holds $N_s \sim \text{Poisson}(\lambda_s)$ animals; each is detected and assigned to a distance bin with probability $\pi_b$, or missed with probability $\pi_0 = 1 - \sum_b \pi_b$. You observe the bin counts $y_{s1}, \ldots, y_{sB}$ and never $N_s$.
Writing the likelihood the obvious way means summing over every possible $N_s$ from the number you saw up to infinity. Royle (2004) noticed you do not have to. Poisson thinning says the bin counts are independent Poisson variables:
$$y_{sb} \sim \text{Poisson}(\lambda_s \pi_b) \quad \text{independently across } b$$
No latent variable, no sum, no truncation ceiling. That is a strong enough claim to be worth checking by brute force rather than by algebra.
```{r identity}
lam <- 4.7; s <- 0.043; y <- c(3, 2, 1, 1, 0)
n <- sum(y); P <- pi_b(s); p0 <- 1 - sum(P); Nmax <- n + 400
closed <- sum(dpois(y, lam*P, log = TRUE))
summed <- log(sum(sapply(n:Nmax, function(N)
dpois(N, lam) *
exp(lfactorial(N) - sum(lfactorial(y)) - lfactorial(N - n)) * # multinomial coefficient
prod(P^y) * p0^(N - n))))
sprintf("latent-N sum to %d: %.12f | independent Poissons: %.12f | difference %.1e",
Nmax, summed, closed, summed - closed)
```
Agreement to `r sprintf("%.0e", abs(summed - closed))`, which is machine arithmetic rather than a good approximation. The identity is exact, and it is exact because the abundance distribution is Poisson. Swap in a negative binomial and it breaks:
```{r counter-world}
th <- 1.5
summed_nb <- log(sum(sapply(n:Nmax, function(N)
dnbinom(N, mu = lam, size = th) *
exp(lfactorial(N) - sum(lfactorial(y)) - lfactorial(N - n)) *
prod(P^y) * p0^(N - n))))
sprintf("N ~ NB(mu = %g, size = %g): true log-likelihood %.6f vs Poisson shortcut %.6f (out by %.4f)",
lam, th, summed_nb, closed, summed_nb - closed)
```
Which is a check worth having in the file: it can fail, and it does fail when it should.
## The world
Density and detectability both respond to the same covariate, which is the ordinary situation rather than a contrived one. Thicker understorey holds more animals and hides them.
$$\log \lambda_s = \beta_0 + \beta_1 z_s, \qquad \log \sigma_s = \alpha_0 + \alpha_1 z_s$$
```{r world}
b0 <- log(4); b1 <- 0.6; a0 <- log(0.05); a1 <- -0.35
sim_hds <- function(seed, alpha1 = a1) {
set.seed(seed)
z <- rnorm(R); lam <- exp(b0 + b1*z); sg <- exp(a0 + alpha1*z)
N <- rpois(R, lam)
y <- t(sapply(seq_len(R), function(i) {
if (N[i] == 0) return(rep(0, B))
x <- runif(N[i], 0, w); d <- x[runif(N[i]) < g_hn(x, sg[i])]
tabulate(findInterval(d, cut, rightmost.closed = TRUE), B)
}))
list(y = y, z = z, N = N)
}
nll_hds <- function(p, y, z, sig_cov = TRUE) {
lam <- exp(p[1] + p[2]*z)
sg <- if (sig_cov) exp(p[3] + p[4]*z) else rep(exp(p[3]), length(z))
-sum(dpois(y, lam * t(sapply(sg, pi_b)), log = TRUE)) # y is R x B, so t() the R x B back
}
fit <- function(par, ...) {
o <- optim(par, nll_hds, ..., method = "Nelder-Mead")
o <- optim(o$par, nll_hds, ..., method = "Nelder-Mead", hessian = TRUE)
list(par = o$par, nll = o$value,
se = tryCatch(sqrt(diag(solve(o$hessian))), error = function(e) rep(NA_real_, length(par))))
}
d <- sim_hds(4282)
sprintf("seed 4282: %d animals present, %d detected (%.1f%%), %d of %d sites detected nothing",
sum(d$N), sum(d$y), 100*sum(d$y)/sum(d$N), sum(rowSums(d$y) == 0), R)
```
At `r sprintf("%.0f", 100*sum(d$y)/sum(d$N))` per cent detected, this is not a hard survey by the standards of the literature. Now fit the full model, and fit the version that models abundance carefully and treats detection as a constant.
```{r fits}
f4 <- fit(c(1, 0.5, log(0.05), -0.3), y = d$y, z = d$z)
f3 <- fit(c(1, 0.5, log(0.05)), y = d$y, z = d$z, sig_cov = FALSE)
rbind(FULL = c(b0 = f4$par[1], b1 = f4$par[2], a0 = f4$par[3], a1 = f4$par[4]),
truth = c(b0, b1, a0, a1))
sprintf("FULL : b1 %.4f (SE %.4f) | a1 %.4f (SE %.4f) | N total %.0f (true %d, raw count %d)",
f4$par[2], f4$se[2], f4$par[4], f4$se[4],
sum(exp(f4$par[1] + f4$par[2]*d$z)), sum(d$N), sum(d$y))
sprintf("NAIVE: b1 %.4f (SE %.4f, %+.1f%%) | 95%% CI [%.3f, %.3f] contains 0.6? %s | dAIC %.1f",
f3$par[2], f3$se[2], 100*(f3$par[2]/b1 - 1),
f3$par[2] - 1.96*f3$se[2], f3$par[2] + 1.96*f3$se[2],
f3$par[2] - 1.96*f3$se[2] < b1 && b1 < f3$par[2] + 1.96*f3$se[2],
(2*f3$nll + 6) - (2*f4$nll + 8))
```
One visit. Not three, not a mark, not a second observer: one pass with a rangefinder, and both the abundance slope and the detection slope come back with usable standard errors.
The naive model returns a habitat effect `r sprintf("%.0f", abs(100*(f3$par[2]/b1 - 1)))` per cent too small, with a confidence interval that excludes the truth. This is the [pooling trade](../covariates-in-the-detection-function/) in hierarchical clothing: detection falls where density rises, and a model with nowhere to put the detection effect books it as fewer animals.
## The exact version of that
One survey is one draw. The asymptotic value the naive slope converges to can be computed directly, by maximising the expected log-likelihood under the true world over a grid of $z$.
```{r pseudo-true}
zg <- seq(-5, 5, length.out = 601); wz <- dnorm(zg); wz <- wz/sum(wz)
Ey <- t(sapply(seq_along(zg), function(i) exp(b0 + b1*zg[i])*pi_b(exp(a0 + a1*zg[i]))))
negQ <- function(p, sig_cov = FALSE) {
lm_ <- exp(p[1] + p[2]*zg)
sg <- if (sig_cov) exp(p[3] + p[4]*zg) else rep(exp(p[3]), length(zg))
M <- t(sapply(seq_along(zg), function(i) lm_[i]*pi_b(sg[i])))
Q <- Ey*log(M); Q[Ey == 0] <- 0 # 0 log 0 = 0, else the truth scores as undefined
-sum(wz * rowSums(Q - M))
}
tight <- function(par, ...) {
o <- optim(par, negQ, ..., method = "Nelder-Mead",
control = list(reltol = 1e-14, maxit = 5000))
for (i in 1:2) o <- optim(o$par, negQ, ..., method = "Nelder-Mead",
control = list(reltol = 1e-14, maxit = 5000))
o$par
}
o4 <- tight(c(1, 0.5, log(0.05), -0.3), sig_cov = TRUE)
o3 <- tight(c(1, 0.5, log(0.05)))
sprintf("FULL model recovers the truth: largest error %.1e", max(abs(o4 - c(b0, b1, a0, a1))))
sprintf("NAIVE converges to: b1 %.4f (truth %.2f, %+.1f%%) | sigma %.5f (true range %.4f to %.4f)",
o3[2], b1, 100*(o3[2]/b1 - 1), exp(o3[3]), exp(a0 + a1*2), exp(a0 - a1*2))
sprintf("implied mean abundance per site: truth %.4f | naive %.4f (%+.1f%%)",
sum(wz*exp(b0 + b1*zg)), sum(wz*exp(o3[1] + o3[2]*zg)),
100*(sum(wz*exp(o3[1] + o3[2]*zg))/sum(wz*exp(b0 + b1*zg)) - 1))
```
The first line is the check, and it earned its place: the correctly specified model must recover the truth exactly, because the expected log-likelihood is maximised at the truth for every value of $z$ at once. It initially did not, by two per cent, which is how the `0 log 0` in the line above got found. A check that cannot fail is decoration.
With that fixed, the naive slope has an exact target: `r sprintf("%.4f", o3[2])`, or `r sprintf("%.0f", abs(100*(o3[2]/b1 - 1)))` per cent low. Note what is *not* badly wrong: mean abundance is out by only `r sprintf("%.1f", abs(100*(sum(wz*exp(o3[1] + o3[2]*zg))/sum(wz*exp(b0 + b1*zg)) - 1)))` per cent. The total survives, the covariate does not, which is the same shape of damage as pooling a detection function across strata, arriving from a different direction.
```{r mc-coverage}
M <- 200
r <- t(sapply(1:M, function(i) {
dd <- sim_hds(52830 + i)
a <- fit(c(1, 0.5, log(0.05), -0.3), y = dd$y, z = dd$z)
b <- fit(c(1, 0.5, log(0.05)), y = dd$y, z = dd$z, sig_cov = FALSE)
c(a$par[2], a$se[2], b$par[2], b$se[2],
sum(exp(a$par[1] + a$par[2]*dd$z))/sum(dd$N),
sum(exp(b$par[1] + b$par[2]*dd$z))/sum(dd$N))
}))
cvg <- function(e, se) mean(abs(e - b1) < 1.96*se)
sprintf("FULL : b1 %.4f (SD %.4f, mean SE %.4f) | CI coverage %.3f | N_hat/N %.4f",
mean(r[, 1]), sd(r[, 1]), mean(r[, 2]), cvg(r[, 1], r[, 2]), mean(r[, 5]))
sprintf("NAIVE: b1 %.4f (SD %.4f, mean SE %.4f) | CI coverage %.3f | N_hat/N %.4f",
mean(r[, 3]), sd(r[, 3]), mean(r[, 4]), cvg(r[, 3], r[, 4]), mean(r[, 6]))
```
The naive slope lands on its asymptotic target of `r sprintf("%.3f", o3[2])` with a standard deviation of `r sprintf("%.3f", sd(r[, 3]))`, and a nominal 95 per cent interval covers the truth in `r sprintf("%.0f", 100*cvg(r[, 3], r[, 4]))` out of `r M` surveys. Not a low coverage: zero. The intervals are the right width for a quantity that is not the one in the title.
```{r fig-fit}
#| echo: false
#| fig-width: 8
#| fig-height: 3.6
#| fig-cap: "Left: observed bin counts against the fitted expectations from the full model, split by whether the site sits above or below the covariate mean. Right: two hundred surveys, showing where each model's habitat slope lands."
#| fig-alt: "Left panel shows observed and fitted distance bin counts for low and high covariate sites, with the high sites having a steeper decline. Right panel shows two histograms of the estimated habitat slope, the full model centred on 0.6 and the naive model centred near 0.32."
ink <- "#16241d"; forest <- "#275139"; gold <- "#c9b458"; rust <- "#b5534e"
hi <- d$z > 0
obs <- rbind(data.frame(bin = 1:B, n = colSums(d$y[!hi, ]), k = "low cover", src = "observed"),
data.frame(bin = 1:B, n = colSums(d$y[hi, ]), k = "high cover", src = "observed"))
lam_h <- exp(f4$par[1] + f4$par[2]*d$z); sg_h <- exp(f4$par[3] + f4$par[4]*d$z)
Efit <- t(sapply(seq_len(R), function(i) lam_h[i]*pi_b(sg_h[i])))
fitd <- rbind(data.frame(bin = 1:B, n = colSums(Efit[!hi, ]), k = "low cover", src = "fitted"),
data.frame(bin = 1:B, n = colSums(Efit[hi, ]), k = "high cover", src = "fitted"))
p1 <- ggplot(obs, aes(bin, n, fill = k)) +
geom_col(position = "dodge", alpha = 0.55) +
geom_line(data = fitd, aes(bin, n, colour = k), linewidth = 0.9) +
geom_point(data = fitd, aes(bin, n, colour = k), size = 1.6) +
scale_fill_manual(values = c("low cover" = forest, "high cover" = rust)) +
scale_colour_manual(values = c("low cover" = forest, "high cover" = rust), guide = "none") +
scale_x_continuous(breaks = 1:B, labels = paste0(cut[-6]*1000, "-", cut[-1]*1000)) +
labs(x = "distance bin (m)", y = "detections", fill = NULL,
title = "Bars observed, lines fitted",
subtitle = "one visit to each of 200 segments") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom", plot.title = element_text(colour = ink, face = "bold"),
plot.subtitle = element_text(size = 8), axis.text.x = element_text(size = 6))
p2 <- ggplot(rbind(data.frame(v = r[, 1], k = "full model"),
data.frame(v = r[, 3], k = "detection held constant")),
aes(v, fill = k)) +
geom_histogram(bins = 32, position = "identity", alpha = 0.6) +
geom_vline(xintercept = b1, colour = ink, linetype = 2, linewidth = 0.7) +
scale_fill_manual(values = c("full model" = forest, "detection held constant" = gold)) +
labs(x = "estimated habitat slope on abundance", y = "surveys", fill = NULL,
title = "Where the slope lands", subtitle = "dashed: the truth (0.6)") +
theme_minimal(base_size = 10) +
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))
```
## One visit with distances against three visits of counts
Now the comparison that decides which survey to run. Same sites, same abundance model, detection held constant at $\sigma = 0.03$ so that the two designs are estimating the same thing. The distance sampler visits each site once and measures. The N-mixture surveyor visits three times and counts.
```{r head-to-head-setup}
s0 <- 0.03; Tv <- 3
pbar <- mu_hn(s0, w)/w
Eg2 <- mu_hn(s0/sqrt(2), w)/w # E[g^2] under x ~ U(0, w)
sprintf("per-visit detection probability p = mu/w = %.5f", pbar)
sprintf("but g(x) across individuals: mean %.4f, SD %.4f -> intra-class correlation %.4f",
pbar, sqrt(Eg2 - pbar^2), (Eg2 - pbar^2)/(pbar*(1 - pbar)))
```
That second line is the whole comparison. The N-mixture model assumes every animal at a site is detected with the same probability $p$ on every visit. In a spatial world there is no such thing as *the* detection probability: an animal on the transect line is seen every time, an animal at 90 m is seen never, and the spread across individuals is enormous.
Whether that breaks the N-mixture depends on something the counts cannot tell you. If the animals hold their positions between visits, as territorial birds do, then each animal carries its own fixed detection probability and the counts are far too consistent between visits. If they shuffle, each animal draws a fresh distance every visit, and the marginal count is exactly $\text{Binomial}(N, \bar{p})$ with $\bar{p} = E[g]$: the N-mixture's assumption holds exactly. Both worlds have the same animals, the same detection function and the same expected count.
```{r head-to-head}
sim_world <- function(seed, territorial) {
set.seed(seed)
z <- rnorm(R); N <- rpois(R, exp(b0 + b1*z))
y1 <- matrix(0, R, B); cnt <- matrix(0, R, Tv)
for (i in seq_len(R)) {
if (N[i] == 0) next
x <- runif(N[i], 0, w) # where the animals are
for (t in seq_len(Tv)) {
xt <- if (territorial) x else runif(N[i], 0, w)
det <- runif(N[i]) < g_hn(xt, s0)
cnt[i, t] <- sum(det)
if (t == 1) y1[i, ] <- tabulate(findInterval(xt[det], cut, rightmost.closed = TRUE), B)
}
}
list(y1 = y1, cnt = cnt, z = z, N = N)
}
## N-mixture likelihood: the binomial coefficients do not move during the fit, so build once
nmix_setup <- function(cnt, Kmax = 90) {
Ng <- 0:Kmax
list(logC = sapply(Ng, function(N) rowSums(lchoose(N, cnt))),
n = rowSums(cnt), Tv = ncol(cnt),
Nmat = matrix(Ng, nrow(cnt), Kmax + 1, byrow = TRUE))
}
nll_nmix <- function(p, S, z) {
lam <- exp(p[1] + p[2]*z); pp <- plogis(p[3])
LP <- dpois(S$Nmat, lam, log = TRUE) + S$logC +
S$n*log(pp) + (S$Tv*S$Nmat - S$n)*log1p(-pp)
-sum(log(rowSums(exp(LP))))
}
run <- function(seed, territorial) {
dw <- sim_world(seed, territorial)
fh <- fit(c(1, 0.5, log(0.04)), y = dw$y1, z = dw$z, sig_cov = FALSE)
S <- nmix_setup(dw$cnt)
fn <- optim(c(1, 0.5, qlogis(0.4)), nll_nmix, S = S, z = dw$z, method = "Nelder-Mead")
fn <- optim(fn$par, nll_nmix, S = S, z = dw$z, method = "Nelder-Mead")
c(hds_b1 = fh$par[2], hds_N = sum(exp(fh$par[1] + fh$par[2]*dw$z))/sum(dw$N),
nm_b1 = fn$par[2], nm_p = plogis(fn$par[3]),
nm_N = sum(exp(fn$par[1] + fn$par[2]*dw$z))/sum(dw$N),
n1 = sum(dw$y1), n3 = sum(dw$cnt),
vv = mean(apply(dw$cnt, 1, var)))
}
for (terr in c(TRUE, FALSE)) {
a <- run(4282, terr)
cat(sprintf("\n--- %s world (seed 4282): %.0f detections in 1 visit, %.0f across 3 visits\n",
if (terr) "TERRITORIAL: animals hold position" else "MOBILE: animals redraw position",
a["n1"], a["n3"]))
cat(sprintf(" distances, 1 visit : b1 %.4f | N_hat/N %.4f\n", a["hds_b1"], a["hds_N"]))
cat(sprintf(" counts, 3 visits : b1 %.4f | p_hat %.4f (true %.4f) | N_hat/N %.4f\n",
a["nm_b1"], a["nm_p"], pbar, a["nm_N"]))
}
```
```{r head-to-head-mc}
M2 <- 100
mc <- lapply(c(TRUE, FALSE), function(terr)
t(sapply(1:M2, function(i) run(52840 + i + 1000*terr, terr))))
names(mc) <- c("territorial", "mobile")
for (k in names(mc)) {
x <- mc[[k]]
cat(sprintf("\n%s world, %d surveys: %.0f detections in 1 visit vs %.0f across 3\n",
k, M2, mean(x[, "n1"]), mean(x[, "n3"])))
cat(sprintf(" distances, 1 visit : b1 %.4f (SD %.4f) | N_hat/N %.4f (SD %.4f)\n",
mean(x[, "hds_b1"]), sd(x[, "hds_b1"]), mean(x[, "hds_N"]), sd(x[, "hds_N"])))
cat(sprintf(" counts, 3 visits : b1 %.4f (SD %.4f) | p_hat %.4f | N_hat/N %.4f (SD %.4f)\n",
mean(x[, "nm_b1"]), sd(x[, "nm_b1"]), mean(x[, "nm_p"]),
mean(x[, "nm_N"]), sd(x[, "nm_N"])))
}
sprintf("mean within-site variance across the 3 visits: territorial %.4f | mobile %.4f",
mean(mc$territorial[, "vv"]), mean(mc$mobile[, "vv"]))
sprintf("territorial N_hat/N should be pbar/p_hat = %.4f; simulated %.4f",
pbar/mean(mc$territorial[, "nm_p"]), mean(mc$territorial[, "nm_N"]))
```
Three results, and only the first one is the one people expect.
**The N-mixture loses `r sprintf("%.0f", 100*(1 - mean(mc$territorial[, "nm_N"])))` per cent of the population in the territorial world**, and it loses it upwards: the fitted detection probability comes out at `r sprintf("%.3f", mean(mc$territorial[, "nm_p"]))` against a truth of `r sprintf("%.3f", pbar)`. The mechanism is the within-site variance in that last line. When animals hold still, the same birds are detected every visit and the counts barely move between visits, `r sprintf("%.2f", mean(mc$territorial[, "vv"]))` against `r sprintf("%.2f", mean(mc$mobile[, "vv"]))`. Consistency is what a binomial reads as high detection, and high detection means there is nothing much left to find. The arithmetic is closed: $\hat{N}/N \to \bar{p}/\hat{p}$, which is `r sprintf("%.4f", pbar/mean(mc$territorial[, "nm_p"]))` against a simulated `r sprintf("%.4f", mean(mc$territorial[, "nm_N"]))`.
Note the direction. Unmodelled individual heterogeneity in capture-recapture inflates abundance, and that is the version most ecologists carry around. Here it deflates it, because the heterogeneity is fixed across visits rather than resampled within them, which is what the between-visit variance is measuring.
**Even in the mobile world, where the N-mixture is unbiased, it is `r sprintf("%.1f", sd(mc$mobile[, "nm_N"])/sd(mc$mobile[, "hds_N"]))` times noisier**, with three times the detections. Repeat visits are an expensive way to buy detection, because they measure it indirectly through a variance.
**The habitat slope is fine everywhere**, near `r sprintf("%.2f", mean(mc$mobile[, "nm_b1"]))` in all four cells. A detection bias that is constant across sites moves the intercept and leaves the slope alone. If your paper is about which habitat has more birds, the N-mixture will not embarrass you. If a number of animals appears in the abstract, it might.
```{r fig-h2h}
#| echo: false
#| fig-width: 8
#| fig-height: 3.6
#| fig-cap: "One hundred surveys in each world. Left: estimated abundance as a fraction of the truth. Right: the fitted per-visit detection probability, which is what the N-mixture gets wrong when animals hold still."
#| fig-alt: "Left panel shows abundance ratios: distance sampling sits on one in both worlds, while the N-mixture sits on one in the mobile world and near 0.56 in the territorial world. Right panel shows the fitted detection probability, correct in the mobile world and inflated in the territorial world."
dfN <- do.call(rbind, lapply(names(mc), function(k) rbind(
data.frame(v = mc[[k]][, "hds_N"], world = k, m = "distances, 1 visit"),
data.frame(v = mc[[k]][, "nm_N"], world = k, m = "counts, 3 visits"))))
p1 <- ggplot(dfN, aes(world, v, fill = m)) +
geom_hline(yintercept = 1, colour = ink, linetype = 2) +
geom_boxplot(alpha = 0.6, outlier.size = 0.5, linewidth = 0.35) +
scale_fill_manual(values = c("distances, 1 visit" = forest, "counts, 3 visits" = gold)) +
labs(x = NULL, y = "estimated N / true N", fill = NULL,
title = "Same animals, same detection function",
subtitle = "dashed: no bias") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom", plot.title = element_text(colour = ink, face = "bold"),
plot.subtitle = element_text(size = 8))
dfp <- do.call(rbind, lapply(names(mc), function(k)
data.frame(v = mc[[k]][, "nm_p"], world = k)))
p2 <- ggplot(dfp, aes(v, fill = world)) +
geom_histogram(bins = 30, position = "identity", alpha = 0.6) +
geom_vline(xintercept = pbar, colour = ink, linetype = 2, linewidth = 0.7) +
scale_fill_manual(values = c(territorial = rust, mobile = forest)) +
labs(x = "fitted per-visit detection probability", y = "surveys", fill = NULL,
title = "What the N-mixture thinks p is",
subtitle = sprintf("dashed: the truth (%.3f)", pbar)) +
theme_minimal(base_size = 10) +
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 this does not fix
The comparison above is not a claim that distance sampling is the better method. It is a claim that the two designs bet on different biology, and that the bets are worth stating out loud.
The N-mixture bets that your animals shuffle between visits, or near enough that a common $p$ is a fair summary. For a territorial songbird in June that bet is poor, and Barker et al. (2018) show the same family of failures from several directions. For a foraging flock it may be fine.
Hierarchical distance sampling bets instead that animals are uniformly distributed with respect to the transect, that everything on the line is detected, and that they hold still *within* a visit rather than moving in response to you. Those are different assumptions, not weaker ones. They are also mostly design choices: random placement makes the first one true by construction, which is why it is worth being fussy about.
And the model still needs the Poisson. The identity that makes the likelihood two lines long is exactly the property that fails under a negative binomial, and unmodelled extra-Poisson variation in abundance has to be added back by hand, at which point the latent-N sum returns along with its truncation ceiling. Sillett et al. (2012) work a full hierarchical distance sampling analysis on an island endemic and are explicit about which of these are assumptions and which are protocol.
The [closing post](../checking-a-distance-sampling-model/) asks what a fitted detection function will actually tell you when one of these is wrong, and the answer is worse than you would guess.
## Related tutorials
- [N-mixture models for abundance](../n-mixture-models-abundance/)
- [Covariates in the detection function](../covariates-in-the-detection-function/)
- [Point transect distance sampling](../point-transect-distance-sampling/)
- [Checking a distance sampling model](../checking-a-distance-sampling-model/)
## References
Barker, R. J., Schofield, M. R., Link, W. A. and Sauer, J. R. (2018). On the reliability of N-mixture models for count data. *Biometrics* 74(1), 369-377. DOI: 10.1111/biom.12734
Buckland, S. T., Rexstad, E. A., Marques, T. A. and Oedekoven, C. S. (2015). *Distance Sampling: Methods and Applications*. Springer. ISBN 978-3-319-19218-5
Kery, M. and Royle, J. A. (2016). *Applied Hierarchical Modeling in Ecology*. Academic Press. ISBN 978-0-12-801378-6
Royle, J. A. (2004). N-mixture models for estimating population size from spatially replicated counts. *Biometrics* 60(1), 108-115. DOI: 10.1111/j.0006-341X.2004.00142.x
Royle, J. A., Dawson, D. K. and Bates, S. (2004). Modeling abundance effects in distance sampling. *Ecology* 85(6), 1591-1597.
Sillett, T. S., Chandler, R. B., Royle, J. A., Kery, M. and Morrison, S. A. (2012). Hierarchical distance-sampling models to estimate population size and habitat-specific abundance of an island endemic. *Ecological Applications* 22(7), 1997-2006.
