Absence of Evidence: How Many Surveys Before You Can Say a Species Is Gone

What a blank notebook page is worth

Somebody walks a transect in the valley, looking for a frog that has not been recorded there since 2011. They find nothing. What have they learned?

Less than it feels like. Two different worlds produce the same blank page, and in one of them the frog is gone. In the other it sits under a rock nine metres away while the surveyor walks past, and the notebook cannot separate those two worlds, nor can any amount of care on the day, because the problem is structural rather than personal. A single visit to an occupied site misses the species with probability \(1 - p\), where \(p\) is whatever the chance of detection happens to be for that species, that habitat and that method.

Repetition is the way out. Each extra blank visit is another chance the species had to show itself and did not take, and at a site where the species really lives the chance of an unbroken run of blanks falls off as the visits accumulate. How fast does it fall off? How many visits does that buy you? And what happens to the number when the assumption holding the whole calculation up, that every visit carries the same detection probability, stops being true, as it always does?

Say plainly what we did. The club has no field site and no licence to handle amphibians, nobody in it has surveyed a transect, and everything below is arithmetic plus simulated survey histories from a seeded generator. The computation is the experiment. Where the model simplifies the world we have tried to name the parts it threw away.

The arithmetic

Take a site where the species is present. Each visit detects it with probability \(p\), independently of the other visits. One blank visit has probability \(1-p\). A run of \(k\) blanks has probability

Set that equal to a target rate \(\alpha\) and solve for \(k\):

With \(\alpha = 0.05\), \(\ln \alpha = -2.9957\). At \(p = 0.20\), \(\ln(0.80) = -0.2231\). The quotient is 13.43. Round up: 14. At \(p = 0.05\) the quotient is 58.40, so 59. At \(p = 0.80\) it is 1.86, so 2.

The ceiling is not cosmetic. Nobody makes 13.43 visits, and 13 of them leave a rate of 0.0550, above the target, which makes rounding up the whole difference between meeting the standard and missing it. Figure 1 shows the decay curves, and Figure 2 shows \(k^{*}\) against \(p\), with the staircase the ceiling produces.

Three properties of the formula matter later. Geometric decay: every visit multiplies the remaining doubt by the same factor. Bare dependence, too, since \(p\) and \(k\) decide the answer between them, with no reference at all to how many sites you hold or how long the season ran. And the condition, which is the property people drop, since the whole expression assumes the site is occupied and therefore says nothing yet about whether the species is absent. Section 10 picks that up.

The model, and the parts of the world it throws away

What we have written down is the detection half of a single-season occupancy model, the framework MacKenzie and colleagues set out in 2002 [1] and that Tyre and colleagues reached from the habitat-modelling side the year after [5]. A small model. Here is what it assumes, and what each assumption quietly deletes.

The site is occupied. Every number in sections 2 through 9 is conditional on that. We compute how often an occupied site produces a blank record, which is a statement about the survey method rather than about the species.

Nothing changes during the survey. Ecologists call this closure. A frog that arrives at the pond halfway through the season, or leaves, breaks it. Real surveys run over weeks, and closure gets violated in ways the data cannot show you.

Visits are independent. We trust this one least. A surveyor who learns the site gets better at finding things, and a disturbed animal hides for a while afterwards. A fortnight of cold spoils four consecutive visits at once, and a wet spring spoils a season, so blanks cluster, and clustered blanks carry less information than the formula credits them with.

Detection probability does not depend on anything. Not the season, not the observer, not the time of day, not the abundance of the animal at that site, and not the number of visits that have already happened there. Royle and Nichols pointed out that abundance alone guarantees the assumption is false, because a site holding twelve individuals is easier to survey than a site holding one, and that pushes detection probability around from site to site whether anyone models it or not [6]. Sections 7 and 8 are about exactly that failure.

No false positives. A detection is real, we assume, though misidentification exists, and for a rare species one wrong identification can do more damage than a hundred missed ones [15].

A model this thin still earns its keep. It gives an exact answer that can be checked, and the places it fails turn out to be the interesting part. Carry one number from this article into a real survey plan, and carry the list above with it.

How we ran it, and what came out

The whole study is one Python file of about 840 lines, using numpy and nothing else, and every random draw descends from a single master seed, 20240921, split into independent streams with numpy's SeedSequence. Each configuration gets its own stream, so rerunning reproduces the printed output exactly.

One replicate is one occupied site given \(k_{\max}\) visits. For each site we draw a row of uniform random numbers, compare each against the detection probability, and record the visit at which the first detection happened, or that none did. The histogram then yields the blank-run probability for every \(k\) at once. We used 500,000 replicates per configuration, ten times the 50,000 minimum we set ourselves, and 2,000,000 for the convergence study, with detection probabilities running from 0.05 to 0.80.

The headline table is below. The "z" column gives the gap between the simulated rate and the closed form, in Monte Carlo standard errors. For a proportion estimated from \(N\) draws that error is \(\sqrt{q(1-q)/N}\). Near zero, the two agree. Above about 3, something is wrong.

Look at the fourth column rather than the third. Rounding \(k^{*}\) up drags the realised false-absence rate well under the 0.05 you asked for, and how far under depends on where the ceiling landed. At \(p = 0.50\) the exact requirement is 4.32 visits, the ceiling gives 5, and the realised rate is 0.0313. You bought far more confidence than you specified. At \(p = 0.05\) the ceiling barely moves anything and the rate lands at 0.0485. Costing a survey, you are looking at free confidence down at low \(p\) and a wasted visit up at high \(p\).

Does the simulator agree with the formula?

Here is what makes the study science rather than output. The simple case has an exact answer, so every simulated number sits beside the number it is supposed to be. We probed nine values of \(k\) at each of eleven detection probabilities, kept the 76 cells where the expected count of blank runs was at least 25, and computed \(z\) for each.

That standard deviation is the line we like best in the whole output, and we tuned nothing to produce it. Two different processes would not hand you a spread of 1.00 by accident.

One subtlety cost us an afternoon: our first summary averaged z over all 76 cells and treated the result as though it carried standard error \(1/\sqrt{76}\). It does not. The nine cells in one row all come from the same 500,000 simulated sites, so a batch that happens to miss slightly too often drags every \(k\) in that row the same way. The cells are strongly correlated, and the pooled mean is far less precise than the count advertises. The fix is one independent anchor cell per detection probability, which gives eleven genuinely independent numbers.

Monte Carlo estimates also have to be shown converging. Otherwise nobody can tell whether a number that agrees agrees for the right reason. Figure 3 tracks the running estimate of \((1-0.20)^{13} = 0.0549756\) as replicates accumulate from 25,000 to 2,000,000, against an envelope of plus or minus two standard errors that narrows as \(1/\sqrt{N}\).

The estimate stays inside the envelope at all fifteen checkpoints, though notice how slowly it improves. Going from 25,000 replicates to 2,000,000, eighty times the work, shrank the standard error from 0.00144 to 0.00016. A factor of nine. That ratio is the square root of eighty, and it is why we care about holding a closed form at all: simulation confirms the formula, and then the formula does the work.

Notes from the club table

Two kinds of variation, and only one of them costs you

Detection probability is not one number. It moves with weather, with season, with whoever is holding the clipboard, and with how many animals are out there. So what does that variation do to the visit requirement? The answer turns entirely on where the variation lives, and the gap is wider than we expected.

Suppose detection probability is drawn from a Beta distribution with mean \(\bar p\). Two quite different arrangements give you that.

Redrawn every visit. Visit \(j\) has its own \(p_j \sim \text{Beta}(a,b)\), drawn fresh, independent of the others. A wet Tuesday followed by a dry Wednesday. Then

because the visits are independent and the miss probability \(1-p\) is linear in \(p\), so the expectation passes straight through the product. The variation vanishes. The spread of the Beta appears nowhere in the answer. A survey whose detection swings between 0.02 and 0.7 from week to week needs exactly as many visits as one sitting at the average every time.

We did not care to take a four-line proof on trust, so we simulated it: four Beta configurations, standard deviations from 0.087 to 0.283, five values of \(k\) each. Largest \(|z|\) across the twenty cells was 2.52, and the four independent anchors averaged \(+0.73\), which is 1.46 standard errors from zero, so no inflation at any spread we tried.

Drawn once and kept. Now the site owns a detection probability, fixed for its life, drawn from the same Beta, because some ponds really are harder than others and stay that way. Every visit to that pond uses that \(p\), so the blank-run probability stops being a product of expectations and becomes instead the expectation of a product:

Because \(x \mapsto (1-x)^k\) is convex, Jensen's inequality puts this at or above \((1-\bar p)^k\), with equality only when the Beta has no spread at all. Now the variation costs you, and the cost compounds with \(k\), because the sites that survive a long run of blanks are disproportionately the hard ones.

We parameterise the Beta by its mean \(\bar p\) and a concentration \(\kappa = a + b\), so that \(a = \kappa \bar p\) and \(b = \kappa(1-\bar p)\), and the standard deviation comes out at \(\sqrt{\bar p (1-\bar p)/(\kappa+1)}\). Small \(\kappa\) means sites differ a lot. Large \(\kappa\) collapses back to the constant case.

The bill for sites that differ

Fix the mean at \(\bar p = 0.20\). The constant-p formula answers 14 visits. Now let the spread grow. Figure 4 gives the blank-run probability on a logarithmic scale, closed form as lines and simulation as rings, for four concentrations plus the constant case.

The constant-p curve is straight on a log scale, because a straight line there is what geometric decay means. None of the heterogeneous curves are straight. They bend upward and flatten, and the flattening is the whole story. For large \(k\) the closed form behaves like

which is a power law rather than a geometric decay, and that one change of shape is the whole of this section. The \(\kappa = 5\) case reads clearest, because there \(a = 1\) exactly and the closed form collapses to \(b/(b+k) = 4/(4+k)\). Twenty blank visits then leave a false-absence probability of exactly \(1/6\). Against that, the naive \(0.80^{20} = 0.0115\). Simulation on 500,000 sites gave 0.167216, and the closed form 0.1666667, which we are content to call agreement.

Set \(4/(4+k) = 0.05\) and you need \(k = 76\) visits, against 14. Here is the same comparison across the grid we ran.

Stare at the right-hand columns. At \(\bar p = 0.60\) with \(\kappa = 2\), twenty blank visits look like a one-in-ninety-million result under the naive formula. They are actually a one-in-forty-five result, an error of a factor of two million. High mean detection is where the naive formula fails worst. The mechanism is plain: a geometric decay with a large \(p\) plunges so fast that any fraction of hard sites, however small, ends up owning the answer.

Figure 5 collects the inflation factors. Read the left-hand edge with care. At \(\kappa = 1\) and \(\bar p = 0.10\) the Beta shape parameter \(a\) is 0.1, the decay exponent is therefore 0.1, and reaching a 5 percent false-absence rate needs \(5.5 \times 10^{12}\) visits. Not a survey design. The model is saying that with that much site-to-site variation the question cannot be settled by visiting more, because a non-negligible fraction of occupied sites are effectively undetectable and always will be.

Three routes to the same number

The heterogeneous case has no obvious right answer to check against, so we built three routes and made them argue: the simulation, the Beta-function closed form via log-gamma, and an independent numerical integration of \(\mathbb{E}[(1-p)^k]\) over the Beta density.

The first integration we wrote was wrong. It split the interval at one half and substituted the endpoint singularity away before applying 240-point Gauss-Legendre, which cures the singularity at one end and manufactures an infinite derivative at the other whenever the exponent exceeds 1. Polynomial quadrature cannot see that. The results disagreed with the closed form by up to 0.86 on a quantity bounded by 1, entirely in the \(\kappa = 100\) rows where \(b\) is large. We replaced it with tanh-sinh quadrature, which handles algebraic endpoint singularities of either sign without any case analysis, and which we set up to integrate the normalising constant too, so that no gamma function appears in that route at all.

Three snakes in Switzerland

All of the above is self-consistent, which is a different thing from being right about anything outside itself. The external test we wanted was somebody who had published both a measured detection probability and the survey effort derived from it. Then check whether our formula reproduces their arithmetic.

Marc Kéry did exactly that in 2002 [2]. Over five years he made 645 visits to 87 sites in Switzerland, following three snake species: the asp viper Vipera aspis, the smooth snake Coronella austriaca, and the grass snake Natrix natrix. He estimated per-visit detection probabilities, found they depended heavily on population size, and reported how many blank visits a site would need before absence could be inferred with 95 percent confidence.

We took his published detection probabilities for small populations and fed them straight into our equation, with nothing fitted to his data and nothing adjusted afterwards.

Two exact. The smooth snake falls two visits short, and we would rather print that row than quietly drop it. Inverting his 34 says the detection probability behind it lay between 0.0843 and 0.0868, whose middle, 0.0855, rounds to the 0.09 in his table. So the gap looks like rounding, not a disagreement about the model. We cannot prove that from a printed table. We are not going to pretend otherwise.

Kéry's other finding should worry anybody quoting our section 2 number in anger. His detection probabilities were not constant. They climbed with population size, from 0.23 to 0.70 for the asp viper and from 0.09 to 0.56 for the smooth snake. Persistent site-level heterogeneity, measured in the field, in exactly the form section 8 prices as expensive. Royle and Nichols had already shown why it happens [6]. Bailey and colleagues found the same structure in salamanders in the Great Smoky Mountains [12]. Garrard and colleagues built a whole survey protocol around estimating detection rates for plants rather than assuming them [7].

The strongest objection we can make to this

Here is the best argument against taking these numbers seriously, put as strongly as we can put it.

You picked the Beta, and nothing in the data says that was right. The inflation factors in section 8 are not properties of heterogeneity as such; they are properties of Beta heterogeneity at one chosen concentration. Change the shape and the answer changes, sometimes by orders of magnitude, and Link showed in 2003 that the shape of a detection-heterogeneity distribution is essentially unidentifiable from capture data, since different mixing distributions can predict nearly identical observations while implying wildly different conclusions [9]. Quoting 76 visits at \(\kappa = 5\) then carries the same epistemic standing as quoting 76 visits because somebody liked the number.

The objection is largely correct and we have no clean answer. Three partial responses, offered as partial.

First, the direction survives even where the magnitude does not. Jensen's inequality does not care what the mixing distribution is. Any spread in a persistent per-site detection probability makes a run of blanks weaker evidence than the constant-p formula claims. So the sign of the error in every published \(k^{*}\) computed from a mean detection probability is known in advance, and it points the dangerous way, which we enjoy rather more than we should.

Second, the asymptotic behaviour is driven by the left tail near \(p = 0\), the part of the distribution a survey can least afford to guess at and also the part every plausible mixing family disagrees about. Distrust the precise value \(5.5 \times 10^{12}\). Take the qualitative claim seriously anyway: the requirement can become unreachable. A two-point mixture, with some fraction of sites near-undetectable, does worse than the Beta. We did not run that case. Good follow-up for somebody.

Third, and weakest: the Beta is what the literature uses, so our numbers are at least comparable with everybody else's [10], which is a reason to publish them alongside. Not a reason to believe them.

A second objection: the number you computed is not the number anyone wants. Also correct. Everything in sections 2 through 9 is conditional on the site being occupied, while a manager wants \(P(\text{absent} \mid k \text{ blanks})\), and that quantity needs a prior. With prior occupancy probability \(\psi\),

We computed the visits needed to push that posterior below 0.05, and at \(p = 0.20\) it takes 24 visits from \(\psi = 0.9\), 14 from \(\psi = 0.5\), and 4 from \(\psi = 0.1\). At \(p = 0.05\): 101, 58 and 15. The prior does most of the work, which is uncomfortable and true, and which is why the decision-theoretic treatments of this problem [8][13][14] put the survey inside a larger question about what you would do with the answer rather than treating the survey as the whole problem.

At \(\psi = 0.5\) the Bayesian answer nearly coincides with the conditional one, 14 against 14 and 58 against 59, a coincidence of the arithmetic at \(\alpha = 0.05\) rather than a deep fact. Also, we suspect, why the conditional number gets quoted as the posterior one.

A third objection, smaller: you treated \(p\) as known. In practice \(p\) is estimated from the same small dataset, often from very few detections, and its uncertainty propagates into \(k^{*}\) nonlinearly, so that because \(k^{*}\) is convex in \(p\) over the range that matters, averaging \(k^{*}\) over an uncertain \(p\) gives a larger answer than plugging in the mean \(p\). Guillera-Arroita and colleagues treat the estimation problem directly [11]. We did not. Read our numbers as conditional on knowing \(p\) exactly, which nobody ever does.

Reproducing this

Clone or download the repository, then from the project root:

python analysis/finding-rare-species.py >
    analysis/finding-rare-species-output.txt

Requirements are Python 3.12 and numpy; we ran Python 3.12.3 and numpy 2.4.2. Expect about 26 seconds on an ordinary laptop. The script reports its own runtime on the last line, and peak memory stays under about 100 MB, because the simulation is chunked at 25,000 sites at a time.

The output runs to 788 lines and is checked into the repository as finding-rare-species-output.txt. Diff your run against ours. Every stream descends from the master seed 20240921 through SeedSequence, so a matching numpy version should reproduce that file byte for byte, apart from the final timing line. Numbers differing in the last digits while the z-columns still sit near zero mean a different numpy bit generator doing its job, not a fault in the method.

Every figure in this article was drawn from that output file. The interactive model runs the same two models in the browser, visit by visit, for anybody who would rather push the parameters around than read a table.

References