Point transect distance sampling

distance sampling

imperfect detection

abundance

ecology tutorial

Points are not short lines. The estimator reads the slope of your distance distribution at zero, which is why the same wrong model costs more here.

Author

Tidy Ecology

Published

2026-08-26

Most bird surveys are not walked, they are stood in. The observer arrives at a point, waits, records every bird detected within some truncation radius, notes how far away each one was, and moves on. Point transects, not line transects.

The machinery looks identical to the line transect case: a detection function, an integral, a count divided by an area. It is not identical, and the difference is not a detail of the arithmetic. On a line, the estimator reads the height of your distance distribution at zero. At a point, it reads the slope. That single change is why point transects need more thought about the key function, and why the standard advice about model choice is stricter here.

library(ggplot2)

g_hn  <- function(r, s) exp(-r^2/(2*s^2))
mu_hn <- function(s, w) sqrt(2*pi)*s*(pnorm(w, 0, s) - 0.5)   # line:  integral of g
nu_hn <- function(s, w) s^2*(1 - g_hn(w, s))                  # point: integral of r g(r)
g_hr  <- function(r, s, b) 1 - exp(-(r/s)^(-b))
mu_hr <- function(s, b, w) integrate(g_hr, 0, w, s = s, b = b)$value
nu_hr <- function(s, b, w) integrate(function(r) r*g_hr(r, s, b), 0, w)$value

Both integrals have closed forms for the half-normal, which is worth checking rather than believing.

s <- 0.037; w <- 0.1
c(mu_closed = mu_hn(s, w), mu_numeric = integrate(g_hn, 0, w, s = s)$value,
  nu_closed = nu_hn(s, w), nu_numeric = integrate(function(r) r*g_hn(r, s), 0, w)$value)

 mu_closed mu_numeric  nu_closed nu_numeric 
0.04605368 0.04605368 0.00133350 0.00133350

The geometry does half the work

On a line, an animal is at perpendicular distance $x$, and distances are uniform over $[0, w]$ by design. At a point, the animals are spread over a disc, and there is more disc further out: the area of the annulus at radius $r$ grows in proportion to $r$. Available distances are not uniform, they are triangular.

So the density of the distances you actually record picks up a factor of $r$:

\[f(r) = \frac{r\,g(r)}{\nu}, \qquad \nu = \int_0^w r\,g(r)\,\mathrm{d}r\]

against the line transect’s $f(x) = g(x)/\mu$. Two consequences follow immediately, and both of them surprise people looking at their first point transect histogram.

w <- 0.1; K <- 200; D <- 100; s0 <- 0.04
nu0 <- nu_hn(s0, w)
Pa  <- 2*nu0/w^2                 # probability an animal in the circle is detected
EDR <- sqrt(2*nu0)               # effective detection radius
En  <- D*K*2*pi*nu0
fpt <- function(r, s, w) r*g_hn(r, s)/nu_hn(s, w)
sprintf("nu %.7f | P_detect %.4f |  EDR %.5f km | E[N in circles] %.1f | E[n detected] %.1f",
        nu0, Pa, EDR, D*K*pi*w^2, En)

[1] "nu 0.0015297 | P_detect 0.3059 |  EDR 0.05531 km | E[N in circles] 628.3 | E[n detected] 192.2"

sprintf("f(0) = %.4f | mode of f at r = %.5f (sigma = %.3f)",
        fpt(0, s0, w), optimise(function(r) -fpt(r, s0, w), c(0, w))$minimum, s0)

[1] "f(0) = 0.0000 | mode of f at r = 0.04001 (sigma = 0.040)"

First, the histogram of detection distances rises from zero, peaks, and falls. The peak sits exactly at $\sigma$ for a half-normal, and none of that is behaviour. The birds are not avoiding the observer. There is simply almost no area at $r = 0$, so almost no birds to detect there, and $f(0) = 0$ exactly whatever the detection function does.

Second, the natural summary of detectability is no longer a width but a radius. The effective detection radius solves $\pi \cdot \mathrm{EDR}^2 = \pi w^2 \cdot P_a$, giving $\mathrm{EDR} = \sqrt{2\nu}$: the radius of the circle within which you would detect everything, and outside which nothing, to give the same expected count. And the density estimator is

\[\hat{D} = \frac{n}{2\pi K \nu} = \frac{n}{K \cdot \pi \cdot \mathrm{EDR}^2}\]

Two ways to get it wrong

The first is to ignore detection entirely and treat the circle as a plot you have censused.

sprintf("plot count D = n/(K pi w^2): %.2f (%+.1f%%)", D*Pa, 100*(Pa - 1))

[1] "plot count D = n/(K pi w^2): 30.59 (-69.4%)"

That is the same mistake as the strip transect, and nobody defends it. The second is subtler and gets made in practice: fit the detection function to the recorded distances as though they were line transect distances, forgetting the $r$ weighting. The likelihood then treats the shortage of birds near zero as evidence that detection is poor near zero, which it reads as a very wide detection function.

## exact: what does the half-normal converge to if the r weighting is dropped?
m2_line  <- function(s, w) s^2*(1 - w*g_hn(w, s)/mu_hn(s, w))
m2_point <- function(s, w) 2*s^4*(1 - (1 + w^2/(2*s^2))*g_hn(w, s))/nu_hn(s, w)
m2_true  <- m2_point(s0, w)
s_naive  <- uniroot(function(s) m2_line(s, w) - m2_true, c(0.005, 0.5), tol = 1e-12)$root
sprintf("E[r^2] = %.7f -> sigma_hat %.5f (true %.3f, %+.1f%%) -> D %.2f (%+.1f%%)",
        m2_true, s_naive, s0, 100*(s_naive/s0 - 1),
        D*nu0/nu_hn(s_naive, w), 100*(nu0/nu_hn(s_naive, w) - 1))

[1] "E[r^2] = 0.0027404 -> sigma_hat 0.08333 (true 0.040, +108.3%) -> D 42.92 (-57.1%)"

## the untruncated version of the same mistake is exactly a factor of two
sprintf("untruncated: E[r^2]/(2 sigma^2) = %.4f, so sigma_hat = sigma*sqrt(2), nu_hat = 2 nu, D halves",
        m2_point(s0, 10)/(2*s0^2))

[1] "untruncated: E[r^2]/(2 sigma^2) = 1.0000, so sigma_hat = sigma*sqrt(2), nu_hat = 2 nu, D halves"

The scale comes out 108 per cent too large and the density comes out 57 per cent too small. Without truncation the arithmetic is exact and memorable: $E[r^2] = 2\sigma^2$ under the point geometry, the line transect estimator solves for $\hat{\sigma} = \sigma\sqrt{2}$, so $\hat{\nu} = 2\nu$ and the density is halved. Precisely halved.

What the estimator is actually reading

Here is the fact that organises everything else. Write the density estimator in terms of the fitted probability density of the observed distances rather than in terms of $g$.

For a line transect, $f(x) = g(x)/\mu$, and since $g(0) = 1$, $f(0) = 1/\mu$, so

\[\hat{D} = \frac{n}{2L\mu} = \frac{n \, f(0)}{2L}\]

For a point transect, $f(r) = r g(r)/\nu$, so $f'(r) = [g(r) + r g'(r)]/\nu$ and $f'(0) = 1/\nu$, so

\[\hat{D} = \frac{n}{2\pi K \nu} = \frac{n \, f'(0)}{2\pi K}\]

h <- 1e-6
fl <- function(x, s, w) g_hn(x, s)/mu_hn(s, w)
sprintf("LINE  f(0) = %.6f vs 1/mu = %.6f (difference %.1e)",
        fl(0, s0, w), 1/mu_hn(s0, w), fl(0, s0, w) - 1/mu_hn(s0, w))

[1] "LINE  f(0) = 20.197959 vs 1/mu = 20.197959 (difference 0.0e+00)"

sprintf("POINT f'(0) = %.6f vs 1/nu = %.6f (difference %.1e, numeric derivative)",
        (fpt(h, s0, w) - fpt(0, s0, w))/h, 1/nu0,
        (fpt(h, s0, w) - fpt(0, s0, w))/h - 1/nu0)

[1] "POINT f'(0) = 653.722565 vs 1/nu = 653.722565 (difference -2.0e-07, numeric derivative)"

Both estimators are one number read off the fitted curve at the origin, and the whole apparatus of key functions and adjustment terms exists to estimate that one number. On a line it is a height. At a point it is a slope.

Estimating a density at a boundary point is hard. Estimating its derivative at a boundary point is harder, and it is hardest exactly where you have fewest observations, because $f(0) = 0$ means the data thin out precisely where the estimator is looking.

set.seed(4281)
N <- rpois(1, D*K*pi*w^2)
r_all <- sqrt(runif(N))*w                    # uniform on the disc: r = w*sqrt(U)
r <- r_all[runif(N) < g_hn(r_all, s0)]
sprintf("n = %d | inside 0.5 sigma: %d (%.1f%%) | inside sigma: %d (%.1f%%)",
        length(r), sum(r < 0.5*s0), 100*mean(r < 0.5*s0), sum(r < s0), 100*mean(r < s0))

[1] "n = 183 | inside 0.5 sigma: 25 (13.7%) | inside sigma: 78 (42.6%)"

One survey, done properly

ld_pt <- function(r, s, w) log(r) - r^2/(2*s^2) - log(nu_hn(s, w))
ld_ln <- function(x, s, w) -x^2/(2*s^2) - log(mu_hn(s, w))
fit_pt <- function(r, w) exp(optimise(function(p) -sum(ld_pt(r, exp(p), w)),
                                      c(-8, 1), tol = 1e-10)$minimum)
fit_ln <- function(x, w) exp(optimise(function(p) -sum(ld_ln(x, exp(p), w)),
                                      c(-8, 1), tol = 1e-10)$minimum)
s_hat  <- fit_pt(r, w); nu_hat <- nu_hn(s_hat, w)
s_wrong <- fit_ln(r, w)
sprintf("correct: sigma %.5f (true %.3f) -> EDR %.5f -> D %.2f (truth %g)",
        s_hat, s0, sqrt(2*nu_hat), length(r)/(2*pi*K*nu_hat), D)

[1] "correct: sigma 0.03905 (true 0.040) -> EDR 0.05418 -> D 99.23 (truth 100)"

sprintf("r weighting dropped: sigma %.5f -> D %.2f (%+.1f%%) | plot count %.2f (%+.1f%%)",
        s_wrong, length(r)/(2*pi*K*nu_hn(s_wrong, w)),
        100*(length(r)/(2*pi*K*nu_hn(s_wrong, w))/D - 1),
        length(r)/(K*pi*w^2), 100*(length(r)/(K*pi*w^2)/D - 1))

[1] "r weighting dropped: sigma 0.07761 -> D 42.87 (-57.1%) | plot count 29.13 (-70.9%)"

Left panel shows a decreasing half-normal detection curve alongside a hump-shaped density that starts at zero and peaks at sigma. Right panel shows a histogram of detection distances rising then falling, with a well-fitting curve and a badly flat curve from the analysis that ignores the r weighting. — Figure 1: Left: the detection function and the distribution of recorded distances it implies, once the growing area of each annulus is accounted for. Right: one survey of 200 points, with the correct fit and the fit that ignores the r weighting.

The same wrong model, two geometries

Now put the slope-versus-height claim to work. Take a world where the true detection function has a shoulder, which real detection functions usually do: everything close is seen, then detection falls away. A hazard-rate key produces that shape; a half-normal cannot. Fit the half-normal anyway, and compare what it costs on a line against what it costs at a point.

s_hr <- 0.05; b_hr <- 3
mu_t <- mu_hr(s_hr, b_hr, w); nu_t <- nu_hr(s_hr, b_hr, w)
sprintf("hazard-rate truth: g(0.5 sigma) = %.4f, g(sigma) = %.4f | mu %.5f | nu %.7f",
        g_hr(0.5*s_hr, s_hr, b_hr), g_hr(s_hr, s_hr, b_hr), mu_t, nu_t)

[1] "hazard-rate truth: g(0.5 sigma) = 0.9997, g(sigma) = 0.6321 | mu 0.06161 | nu 0.0021177"

best_fit <- function(dens, ld, w) {
  optimise(function(p) -integrate(function(x) dens(x)*ld(x, exp(p), w), 1e-9, w)$value,
           c(-6, 0), tol = 1e-12)$minimum
}
sL <- exp(best_fit(function(x) g_hr(x, s_hr, b_hr)/mu_t, ld_ln, w))
sP <- exp(best_fit(function(r) r*g_hr(r, s_hr, b_hr)/nu_t, ld_pt, w))
sprintf("LINE : best half-normal sigma %.5f -> mu %.5f vs true %.5f -> D bias %+.2f%%",
        sL, mu_hn(sL, w), mu_t, 100*(mu_t/mu_hn(sL, w) - 1))

[1] "LINE : best half-normal sigma 0.04544 -> mu 0.05537 vs true 0.06161 -> D bias +11.27%"

sprintf("POINT: best half-normal sigma %.5f -> nu %.7f vs true %.7f -> D bias %+.2f%%",
        sP, nu_hn(sP, w), nu_t, 100*(nu_t/nu_hn(sP, w) - 1))

[1] "POINT: best half-normal sigma 0.04445 -> nu 0.0018185 vs true 0.0021177 -> D bias +16.45%"

The same key function, fitted to the same detection process, costs 11.3 per cent on a line and 16.5 per cent at a point: about 1.5 times worse. Both biases are positive, because a model with no shoulder fitted to shouldered data underestimates the effective width or radius, and dividing by too small an area inflates density.

The bias is not an approximation of anything. It is the error in that one number at the origin, restated: the ratio $\hat{D}/D$ equals the ratio of fitted to true slope at zero, exactly, and the point geometry makes that slope harder to pin down.

Truncating harder is the usual reflex, and it does not rescue this.

sw <- t(sapply(c(0.06, 0.07, 0.08, 0.09, 0.10), function(ww) {
  nt <- nu_hr(s_hr, b_hr, ww); mt <- mu_hr(s_hr, b_hr, ww)
  sp <- exp(best_fit(function(r) r*g_hr(r, s_hr, b_hr)/nt, ld_pt, ww))
  sl <- exp(best_fit(function(x) g_hr(x, s_hr, b_hr)/mt, ld_ln, ww))
  c(w = ww, dropped_pct = 100*(1 - nt/nu_t),
    point_bias = 100*(nt/nu_hn(sp, ww) - 1), line_bias = 100*(mt/mu_hn(sl, ww) - 1))
}))
round(sw, 2)

        w dropped_pct point_bias line_bias
[1,] 0.06       33.85      16.50      8.14
[2,] 0.07       22.61      20.00     10.45
[3,] 0.08       13.51      20.70     11.58
[4,] 0.09        6.10      19.27     11.77
[5,] 0.10        0.00      16.45     11.27

Cutting the truncation radius back to 60 m throws away a third of the detections and leaves the point transect bias exactly where it started. The bias is not living in the tail.

Four hundred surveys of each

To compare the two designs fairly, give the line survey enough transect to expect the same number of detections as the 200 points, and let both analysts make the same modelling mistake.

Ln <- 266.1/(2*D*mu_t)
M <- 400
res <- t(sapply(1:M, function(i) {
  set.seed(52820 + i)
  Np <- rpois(1, D*K*pi*w^2); rp <- sqrt(runif(Np))*w
  rp <- rp[runif(Np) < g_hr(rp, s_hr, b_hr)]
  Nl <- rpois(1, D*2*w*Ln);   xl <- runif(Nl, 0, w)
  xl <- xl[runif(Nl) < g_hr(xl, s_hr, b_hr)]
  c(n_pt = length(rp), n_ln = length(xl),
    D_pt = length(rp)/(2*pi*K*nu_hn(fit_pt(rp, w), w)),
    D_ln = length(xl)/(2*Ln*mu_hn(fit_ln(xl, w), w)))
}))
sprintf("effort matched: %d points and %.1f km of line, both expecting %.0f detections",
        K, Ln, 2*D*Ln*mu_t)

[1] "effort matched: 200 points and 21.6 km of line, both expecting 266 detections"

sprintf("POINT: D %.2f (SD %.2f) = %+.1f%% | within 5%% of the truth in %.1f%% of surveys",
        mean(res[, "D_pt"]), sd(res[, "D_pt"]), 100*(mean(res[, "D_pt"])/D - 1),
        100*mean(abs(res[, "D_pt"]/D - 1) < 0.05))

[1] "POINT: D 117.81 (SD 10.83) = +17.8% | within 5% of the truth in 11.5% of surveys"

sprintf("LINE : D %.2f (SD %.2f) = %+.1f%% | within 5%% of the truth in %.1f%% of surveys",
        mean(res[, "D_ln"]), sd(res[, "D_ln"]), 100*(mean(res[, "D_ln"])/D - 1),
        100*mean(abs(res[, "D_ln"]/D - 1) < 0.05))

[1] "LINE : D 111.59 (SD 8.61) = +11.6% | within 5% of the truth in 21.8% of surveys"

The simulated biases run a little above the exact asymptotic values (+16.5 and +11.3 per cent), which is what happens when you average a ratio: $E[1/\hat{\nu}]$ sits above $1/E[\hat{\nu}]$. The ordering is the point, and it is stable: equal detections, equal mistake, and the point survey lands near the truth half as often.

Left panel shows a shouldered hazard-rate detection curve with two half-normal curves that miss the shoulder. Right panel shows two histograms of density estimates, both above the true value, with the point transect one further out and wider. — Figure 2: Left: the true hazard-rate curve and the best half-normal approximations under the two geometries, which differ mainly in their slope at the origin. Right: four hundred effort-matched surveys under each design, both analysed with the wrong key function.

What this means in the field

The practical advice that falls out of the slope has been standard since Buckland (2006), and now it has a reason attached rather than a rule to memorise.

Prefer a key function with a shoulder, and be readier to add adjustment terms at a point than on a line. Group the distances rather than pretending 3 m of precision at 90 m, since binned data fitted with a multinomial likelihood are honest about what you measured. Record enough detections: the usual line transect target of 60 to 80 is thin for a point survey. And take the near-point measurements seriously, because that handful of birds inside half a sigma is carrying the slope.

What this does not fix

The estimator assumes the animals are uniformly distributed with respect to the point, which is a statement about where you put your points, not about the birds. Points placed at path junctions or in clearings break it, and no amount of care with the key function repairs a density gradient around the point.

It also assumes everything at the point is detected. On a line an observer walks along and can be reasonably sure of the animals on the line itself. Standing at a point, the birds directly overhead in the canopy are not obviously easier. That assumption is the subject of the checking post, and the short version is that your distances cannot test it.

Finally, this post held detection to a single key function with a single scale. Real point counts vary by observer, by time of morning, by wind and by species, all of which belong in the detection function rather than in the residual.

References

Buckland, S. T. (2006). Point transect surveys for songbirds: robust methodologies. The Auk 123(2), 345-357. DOI: 10.1093/auk/123.2.345

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

Marques, F. F. C. and Buckland, S. T. (2003). Incorporating covariates into standard line transect analyses. Biometrics 59(4), 924-935. DOI: 10.1111/j.0006-341X.2003.00107.x

Marques, T. A., Thomas, L., Fancy, S. G. and Buckland, S. T. (2007). Improving estimates of bird density using multiple-covariate distance sampling. The Auk 124(4), 1229-1243.

Rexstad, E., Buckland, S. T., Marshall, L. and Borchers, D. (2023). Pooling robustness in distance sampling: avoiding bias when there is unmodelled heterogeneity. Ecology and Evolution 13(1), e9684. DOI: 10.1002/ece3.9684

--- title: "Point transect distance sampling" description: "Points are not short lines. The estimator reads the slope of your distance distribution at zero, which is why the same wrong model costs more here." date: 2026-08-26 10:00 categories: [R, distance sampling, imperfect detection, abundance, ecology tutorial] image: thumbnail.png image-alt: "A detection curve beside the distribution of detection distances rising from zero to a peak" --- Most bird surveys are not walked, they are stood in. The observer arrives at a point, waits, records every bird detected within some truncation radius, notes how far away each one was, and moves on. Point transects, not line transects. The machinery looks identical to [the line transect case](../distance-sampling-line-transect/): a detection function, an integral, a count divided by an area. It is not identical, and the difference is not a detail of the arithmetic. On a line, the estimator reads the height of your distance distribution at zero. At a point, it reads the slope. That single change is why point transects need more thought about the key function, and why the standard advice about model choice is stricter here. ```{r setup} #| message: false #| warning: false library(ggplot2) g_hn <- function(r, s) exp(-r^2/(2*s^2)) mu_hn <- function(s, w) sqrt(2*pi)*s*(pnorm(w, 0, s) - 0.5) # line: integral of g nu_hn <- function(s, w) s^2*(1 - g_hn(w, s)) # point: integral of r g(r) g_hr <- function(r, s, b) 1 - exp(-(r/s)^(-b)) mu_hr <- function(s, b, w) integrate(g_hr, 0, w, s = s, b = b)$value nu_hr <- function(s, b, w) integrate(function(r) r*g_hr(r, s, b), 0, w)$value ``` Both integrals have closed forms for the half-normal, which is worth checking rather than believing. ```{r check-closed} s <- 0.037; w <- 0.1 c(mu_closed = mu_hn(s, w), mu_numeric = integrate(g_hn, 0, w, s = s)$value, nu_closed = nu_hn(s, w), nu_numeric = integrate(function(r) r*g_hn(r, s), 0, w)$value) ``` ## The geometry does half the work On a line, an animal is at perpendicular distance $x$, and distances are uniform over $[0, w]$ by design. At a point, the animals are spread over a disc, and there is more disc further out: the area of the annulus at radius $r$ grows in proportion to $r$. Available distances are not uniform, they are triangular. So the density of the distances you actually record picks up a factor of $r$: $$f(r) = \frac{r\,g(r)}{\nu}, \qquad \nu = \int_0^w r\,g(r)\,\mathrm{d}r$$ against the line transect's $f(x) = g(x)/\mu$. Two consequences follow immediately, and both of them surprise people looking at their first point transect histogram. ```{r geometry} w <- 0.1; K <- 200; D <- 100; s0 <- 0.04 nu0 <- nu_hn(s0, w) Pa <- 2*nu0/w^2 # probability an animal in the circle is detected EDR <- sqrt(2*nu0) # effective detection radius En <- D*K*2*pi*nu0 fpt <- function(r, s, w) r*g_hn(r, s)/nu_hn(s, w) sprintf("nu %.7f | P_detect %.4f | EDR %.5f km | E[N in circles] %.1f | E[n detected] %.1f", nu0, Pa, EDR, D*K*pi*w^2, En) sprintf("f(0) = %.4f | mode of f at r = %.5f (sigma = %.3f)", fpt(0, s0, w), optimise(function(r) -fpt(r, s0, w), c(0, w))$minimum, s0) ``` First, the histogram of detection distances rises from zero, peaks, and falls. The peak sits exactly at $\sigma$ for a half-normal, and none of that is behaviour. The birds are not avoiding the observer. There is simply almost no area at $r = 0$, so almost no birds to detect there, and $f(0) = 0$ exactly whatever the detection function does. Second, the natural summary of detectability is no longer a width but a radius. The effective detection radius solves $\pi \cdot \mathrm{EDR}^2 = \pi w^2 \cdot P_a$, giving $\mathrm{EDR} = \sqrt{2\nu}$: the radius of the circle within which you would detect everything, and outside which nothing, to give the same expected count. And the density estimator is $$\hat{D} = \frac{n}{2\pi K \nu} = \frac{n}{K \cdot \pi \cdot \mathrm{EDR}^2}$$ ## Two ways to get it wrong The first is to ignore detection entirely and treat the circle as a plot you have censused. ```{r naive-plot} sprintf("plot count D = n/(K pi w^2): %.2f (%+.1f%%)", D*Pa, 100*(Pa - 1)) ``` That is the same mistake as the strip transect, and nobody defends it. The second is subtler and gets made in practice: fit the detection function to the recorded distances as though they were line transect distances, forgetting the $r$ weighting. The likelihood then treats the shortage of birds near zero as evidence that detection is poor near zero, which it reads as a very wide detection function. ```{r naive-weight} ## exact: what does the half-normal converge to if the r weighting is dropped? m2_line <- function(s, w) s^2*(1 - w*g_hn(w, s)/mu_hn(s, w)) m2_point <- function(s, w) 2*s^4*(1 - (1 + w^2/(2*s^2))*g_hn(w, s))/nu_hn(s, w) m2_true <- m2_point(s0, w) s_naive <- uniroot(function(s) m2_line(s, w) - m2_true, c(0.005, 0.5), tol = 1e-12)$root sprintf("E[r^2] = %.7f -> sigma_hat %.5f (true %.3f, %+.1f%%) -> D %.2f (%+.1f%%)", m2_true, s_naive, s0, 100*(s_naive/s0 - 1), D*nu0/nu_hn(s_naive, w), 100*(nu0/nu_hn(s_naive, w) - 1)) ## the untruncated version of the same mistake is exactly a factor of two sprintf("untruncated: E[r^2]/(2 sigma^2) = %.4f, so sigma_hat = sigma*sqrt(2), nu_hat = 2 nu, D halves", m2_point(s0, 10)/(2*s0^2)) ``` The scale comes out `r sprintf("%.0f", 100*(s_naive/s0 - 1))` per cent too large and the density comes out `r sprintf("%.0f", abs(100*(nu0/nu_hn(s_naive, w) - 1)))` per cent too small. Without truncation the arithmetic is exact and memorable: $E[r^2] = 2\sigma^2$ under the point geometry, the line transect estimator solves for $\hat{\sigma} = \sigma\sqrt{2}$, so $\hat{\nu} = 2\nu$ and the density is halved. Precisely halved. ## What the estimator is actually reading Here is the fact that organises everything else. Write the density estimator in terms of the fitted probability density of the observed distances rather than in terms of $g$. For a line transect, $f(x) = g(x)/\mu$, and since $g(0) = 1$, $f(0) = 1/\mu$, so $$\hat{D} = \frac{n}{2L\mu} = \frac{n \, f(0)}{2L}$$ For a point transect, $f(r) = r g(r)/\nu$, so $f'(r) = [g(r) + r g'(r)]/\nu$ and $f'(0) = 1/\nu$, so $$\hat{D} = \frac{n}{2\pi K \nu} = \frac{n \, f'(0)}{2\pi K}$$ ```{r identities} h <- 1e-6 fl <- function(x, s, w) g_hn(x, s)/mu_hn(s, w) sprintf("LINE f(0) = %.6f vs 1/mu = %.6f (difference %.1e)", fl(0, s0, w), 1/mu_hn(s0, w), fl(0, s0, w) - 1/mu_hn(s0, w)) sprintf("POINT f'(0) = %.6f vs 1/nu = %.6f (difference %.1e, numeric derivative)", (fpt(h, s0, w) - fpt(0, s0, w))/h, 1/nu0, (fpt(h, s0, w) - fpt(0, s0, w))/h - 1/nu0) ``` Both estimators are one number read off the fitted curve at the origin, and the whole apparatus of key functions and adjustment terms exists to estimate that one number. On a line it is a height. At a point it is a slope. Estimating a density at a boundary point is hard. Estimating its derivative at a boundary point is harder, and it is hardest exactly where you have fewest observations, because $f(0) = 0$ means the data thin out precisely where the estimator is looking. ```{r sparse} set.seed(4281) N <- rpois(1, D*K*pi*w^2) r_all <- sqrt(runif(N))*w # uniform on the disc: r = w*sqrt(U) r <- r_all[runif(N) < g_hn(r_all, s0)] sprintf("n = %d | inside 0.5 sigma: %d (%.1f%%) | inside sigma: %d (%.1f%%)", length(r), sum(r < 0.5*s0), 100*mean(r < 0.5*s0), sum(r < s0), 100*mean(r < s0)) ``` ## One survey, done properly ```{r worked} ld_pt <- function(r, s, w) log(r) - r^2/(2*s^2) - log(nu_hn(s, w)) ld_ln <- function(x, s, w) -x^2/(2*s^2) - log(mu_hn(s, w)) fit_pt <- function(r, w) exp(optimise(function(p) -sum(ld_pt(r, exp(p), w)), c(-8, 1), tol = 1e-10)$minimum) fit_ln <- function(x, w) exp(optimise(function(p) -sum(ld_ln(x, exp(p), w)), c(-8, 1), tol = 1e-10)$minimum) s_hat <- fit_pt(r, w); nu_hat <- nu_hn(s_hat, w) s_wrong <- fit_ln(r, w) sprintf("correct: sigma %.5f (true %.3f) -> EDR %.5f -> D %.2f (truth %g)", s_hat, s0, sqrt(2*nu_hat), length(r)/(2*pi*K*nu_hat), D) sprintf("r weighting dropped: sigma %.5f -> D %.2f (%+.1f%%) | plot count %.2f (%+.1f%%)", s_wrong, length(r)/(2*pi*K*nu_hn(s_wrong, w)), 100*(length(r)/(2*pi*K*nu_hn(s_wrong, w))/D - 1), length(r)/(K*pi*w^2), 100*(length(r)/(K*pi*w^2)/D - 1)) ``` ```{r fig-geometry} #| echo: false #| fig-width: 8 #| fig-height: 3.6 #| fig-cap: "Left: the detection function and the distribution of recorded distances it implies, once the growing area of each annulus is accounted for. Right: one survey of 200 points, with the correct fit and the fit that ignores the r weighting." #| fig-alt: "Left panel shows a decreasing half-normal detection curve alongside a hump-shaped density that starts at zero and peaks at sigma. Right panel shows a histogram of detection distances rising then falling, with a well-fitting curve and a badly flat curve from the analysis that ignores the r weighting." ink <- "#16241d"; forest <- "#275139"; gold <- "#c9b458"; rust <- "#b5534e" rs <- seq(0, w, length.out = 400) p1 <- ggplot(data.frame(r = rep(rs, 2), y = c(g_hn(rs, s0), fpt(rs, s0, w)/max(fpt(rs, s0, w))), k = rep(c("g(r): detection probability", "f(r): distances you record"), each = length(rs))), aes(r*1000, y, colour = k)) + geom_line(linewidth = 0.9) + geom_vline(xintercept = s0*1000, colour = ink, linetype = 3) + scale_colour_manual(values = c("g(r): detection probability" = forest, "f(r): distances you record" = gold)) + labs(x = "distance from the point (m)", y = "scaled to a maximum of one", colour = NULL, title = "Area rises, detection falls", subtitle = "dotted: sigma, where f(r) peaks") + theme_minimal(base_size = 10) + theme(legend.position = "bottom", plot.title = element_text(colour = ink, face = "bold"), plot.subtitle = element_text(size = 8)) fitcurves <- rbind( data.frame(r = rs, f = fpt(rs, s_hat, w), k = "correct fit"), data.frame(r = rs, f = fpt(rs, s_wrong, w), k = "r weighting dropped")) p2 <- ggplot(data.frame(r = r), aes(r*1000)) + geom_histogram(aes(y = after_stat(density)/1000), bins = 18, fill = forest, alpha = 0.35, colour = NA) + geom_line(data = fitcurves, aes(r*1000, f/1000, colour = k), linewidth = 0.9) + scale_colour_manual(values = c("correct fit" = forest, "r weighting dropped" = rust)) + labs(x = "distance from the point (m)", y = "density", colour = NULL, title = "One survey, two analyses", subtitle = sprintf("n = %d detections at %d points", length(r), K)) + 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)) ``` ## The same wrong model, two geometries Now put the slope-versus-height claim to work. Take a world where the true detection function has a shoulder, which real detection functions usually do: everything close is seen, then detection falls away. A hazard-rate key produces that shape; a half-normal cannot. Fit the half-normal anyway, and compare what it costs on a line against what it costs at a point. ```{r misspec} s_hr <- 0.05; b_hr <- 3 mu_t <- mu_hr(s_hr, b_hr, w); nu_t <- nu_hr(s_hr, b_hr, w) sprintf("hazard-rate truth: g(0.5 sigma) = %.4f, g(sigma) = %.4f | mu %.5f | nu %.7f", g_hr(0.5*s_hr, s_hr, b_hr), g_hr(s_hr, s_hr, b_hr), mu_t, nu_t) best_fit <- function(dens, ld, w) { optimise(function(p) -integrate(function(x) dens(x)*ld(x, exp(p), w), 1e-9, w)$value, c(-6, 0), tol = 1e-12)$minimum } sL <- exp(best_fit(function(x) g_hr(x, s_hr, b_hr)/mu_t, ld_ln, w)) sP <- exp(best_fit(function(r) r*g_hr(r, s_hr, b_hr)/nu_t, ld_pt, w)) sprintf("LINE : best half-normal sigma %.5f -> mu %.5f vs true %.5f -> D bias %+.2f%%", sL, mu_hn(sL, w), mu_t, 100*(mu_t/mu_hn(sL, w) - 1)) sprintf("POINT: best half-normal sigma %.5f -> nu %.7f vs true %.7f -> D bias %+.2f%%", sP, nu_hn(sP, w), nu_t, 100*(nu_t/nu_hn(sP, w) - 1)) ``` The same key function, fitted to the same detection process, costs `r sprintf("%.1f", 100*(mu_t/mu_hn(sL, w) - 1))` per cent on a line and `r sprintf("%.1f", 100*(nu_t/nu_hn(sP, w) - 1))` per cent at a point: about `r sprintf("%.1f", (nu_t/nu_hn(sP, w) - 1)/(mu_t/mu_hn(sL, w) - 1))` times worse. Both biases are positive, because a model with no shoulder fitted to shouldered data underestimates the effective width or radius, and dividing by too small an area inflates density. The bias is not an approximation of anything. It is the error in that one number at the origin, restated: the ratio $\hat{D}/D$ equals the ratio of fitted to true slope at zero, exactly, and the point geometry makes that slope harder to pin down. Truncating harder is the usual reflex, and it does not rescue this. ```{r truncation-sweep} sw <- t(sapply(c(0.06, 0.07, 0.08, 0.09, 0.10), function(ww) { nt <- nu_hr(s_hr, b_hr, ww); mt <- mu_hr(s_hr, b_hr, ww) sp <- exp(best_fit(function(r) r*g_hr(r, s_hr, b_hr)/nt, ld_pt, ww)) sl <- exp(best_fit(function(x) g_hr(x, s_hr, b_hr)/mt, ld_ln, ww)) c(w = ww, dropped_pct = 100*(1 - nt/nu_t), point_bias = 100*(nt/nu_hn(sp, ww) - 1), line_bias = 100*(mt/mu_hn(sl, ww) - 1)) })) round(sw, 2) ``` Cutting the truncation radius back to 60 m throws away a third of the detections and leaves the point transect bias exactly where it started. The bias is not living in the tail. ## Four hundred surveys of each To compare the two designs fairly, give the line survey enough transect to expect the same number of detections as the 200 points, and let both analysts make the same modelling mistake. ```{r mc} Ln <- 266.1/(2*D*mu_t) M <- 400 res <- t(sapply(1:M, function(i) { set.seed(52820 + i) Np <- rpois(1, D*K*pi*w^2); rp <- sqrt(runif(Np))*w rp <- rp[runif(Np) < g_hr(rp, s_hr, b_hr)] Nl <- rpois(1, D*2*w*Ln); xl <- runif(Nl, 0, w) xl <- xl[runif(Nl) < g_hr(xl, s_hr, b_hr)] c(n_pt = length(rp), n_ln = length(xl), D_pt = length(rp)/(2*pi*K*nu_hn(fit_pt(rp, w), w)), D_ln = length(xl)/(2*Ln*mu_hn(fit_ln(xl, w), w))) })) sprintf("effort matched: %d points and %.1f km of line, both expecting %.0f detections", K, Ln, 2*D*Ln*mu_t) sprintf("POINT: D %.2f (SD %.2f) = %+.1f%% | within 5%% of the truth in %.1f%% of surveys", mean(res[, "D_pt"]), sd(res[, "D_pt"]), 100*(mean(res[, "D_pt"])/D - 1), 100*mean(abs(res[, "D_pt"]/D - 1) < 0.05)) sprintf("LINE : D %.2f (SD %.2f) = %+.1f%% | within 5%% of the truth in %.1f%% of surveys", mean(res[, "D_ln"]), sd(res[, "D_ln"]), 100*(mean(res[, "D_ln"])/D - 1), 100*mean(abs(res[, "D_ln"]/D - 1) < 0.05)) ``` The simulated biases run a little above the exact asymptotic values (`r sprintf("%+.1f", 100*(nu_t/nu_hn(sP, w) - 1))` and `r sprintf("%+.1f", 100*(mu_t/mu_hn(sL, w) - 1))` per cent), which is what happens when you average a ratio: $E[1/\hat{\nu}]$ sits above $1/E[\hat{\nu}]$. The ordering is the point, and it is stable: equal detections, equal mistake, and the point survey lands near the truth half as often. ```{r fig-mc} #| echo: false #| fig-width: 8 #| fig-height: 3.6 #| fig-cap: "Left: the true hazard-rate curve and the best half-normal approximations under the two geometries, which differ mainly in their slope at the origin. Right: four hundred effort-matched surveys under each design, both analysed with the wrong key function." #| fig-alt: "Left panel shows a shouldered hazard-rate detection curve with two half-normal curves that miss the shoulder. Right panel shows two histograms of density estimates, both above the true value, with the point transect one further out and wider." p1 <- ggplot(data.frame(r = rep(rs, 3), g = c(g_hr(rs, s_hr, b_hr), g_hn(rs, sL), g_hn(rs, sP)), k = rep(c("hazard-rate (true)", "best half-normal, line", "best half-normal, point"), each = length(rs))), aes(r*1000, g, colour = k, linetype = k)) + geom_line(linewidth = 0.9) + scale_colour_manual(values = c("hazard-rate (true)" = ink, "best half-normal, line" = forest, "best half-normal, point" = rust)) + scale_linetype_manual(values = c("hazard-rate (true)" = 1, "best half-normal, line" = 2, "best half-normal, point" = 2)) + labs(x = "distance (m)", y = "detection probability", colour = NULL, linetype = NULL, title = "A shoulder the half-normal cannot make") + theme_minimal(base_size = 10) + theme(legend.position = "bottom", plot.title = element_text(colour = ink, face = "bold"), legend.text = element_text(size = 7)) mc <- rbind(data.frame(v = res[, "D_pt"], k = "point (200 points)"), data.frame(v = res[, "D_ln"], k = sprintf("line (%.1f km)", Ln))) p2 <- ggplot(mc, aes(v, fill = k)) + geom_histogram(bins = 34, position = "identity", alpha = 0.55) + geom_vline(xintercept = D, colour = ink, linetype = 2, linewidth = 0.7) + scale_fill_manual(values = setNames(c(rust, forest), c("point (200 points)", sprintf("line (%.1f km)", Ln)))) + labs(x = "estimated density (animals per km2)", y = "surveys", fill = NULL, title = "Equal detections, equal mistake", subtitle = "dashed: the truth (100)") + 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 means in the field The practical advice that falls out of the slope has been standard since Buckland (2006), and now it has a reason attached rather than a rule to memorise. Prefer a key function with a shoulder, and be readier to add adjustment terms at a point than on a line. Group the distances rather than pretending 3 m of precision at 90 m, since binned data fitted with a multinomial likelihood are honest about what you measured. Record enough detections: the usual line transect target of 60 to 80 is thin for a point survey. And take the near-point measurements seriously, because that handful of birds inside half a sigma is carrying the slope. ## What this does not fix The estimator assumes the animals are uniformly distributed with respect to the point, which is a statement about where you put your points, not about the birds. Points placed at path junctions or in clearings break it, and no amount of care with the key function repairs a density gradient around the point. It also assumes everything at the point is detected. On a line an observer walks along and can be reasonably sure of the animals on the line itself. Standing at a point, the birds directly overhead in the canopy are not obviously easier. That assumption is the subject of the [checking post](../checking-a-distance-sampling-model/), and the short version is that your distances cannot test it. Finally, this post held detection to a single key function with a single scale. Real point counts vary by observer, by time of morning, by wind and by species, all of which belong [in the detection function](../covariates-in-the-detection-function/) rather than in the residual. ## Related tutorials - [Distance sampling on line transects](../distance-sampling-line-transect/) - [Covariates in the detection function](../covariates-in-the-detection-function/) - [Hierarchical distance sampling](../hierarchical-distance-sampling/) - [Checking a distance sampling model](../checking-a-distance-sampling-model/) ## References Buckland, S. T. (2006). Point transect surveys for songbirds: robust methodologies. *The Auk* 123(2), 345-357. DOI: 10.1093/auk/123.2.345 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 Marques, F. F. C. and Buckland, S. T. (2003). Incorporating covariates into standard line transect analyses. *Biometrics* 59(4), 924-935. DOI: 10.1111/j.0006-341X.2003.00107.x Marques, T. A., Thomas, L., Fancy, S. G. and Buckland, S. T. (2007). Improving estimates of bird density using multiple-covariate distance sampling. *The Auk* 124(4), 1229-1243. Rexstad, E., Buckland, S. T., Marshall, L. and Borchers, D. (2023). Pooling robustness in distance sampling: avoiding bias when there is unmodelled heterogeneity. *Ecology and Evolution* 13(1), e9684. DOI: 10.1002/ece3.9684