Estimating extinction dates from sighting records

Overview

EDE implements methods for inference about extinction from sighting records and fossil ranges. Some estimate an endpoint directly; others test whether a species could plausibly have persisted to a candidate time. The tests return frequentist p-values, not posterior probabilities that the species is extant; the Bayesian models later in this vignette give those. Their source papers define no point estimate, so EDE reports the one-sided upper confidence bound obtained by inverting each test.

The column containing these p-values is currently named chance for compatibility with EDE 0.1.0.

Input and observation origin

library(EDE)

years <- c(1900, 1902, 1903, 1905, 1907, 1908, 1910, 1912,
           1915, 1918, 1920, 1923, 1925, 1928, 1930, 1933, 1936)
counts <- c(4, 3, 5, 2, 3, 4, 2, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1)
sd <- sighting_data(data.frame(year = years, sightings = counts))

The earliest supplied time defines the observation origin. Include an initial zero-count row if observation began before the first sighting. Methods differ in their treatment of counts, following their source papers: OLE, the Weibull test, both Solow (1993) tests and Burgman’s tests count every sighting event; the Robson-Whitlock, Strauss-Sadler, McInerny, Solow-Roberts and Jarić-Ebenhard methods use only distinct sighting times. Strauss-Sadler, Solow-Roberts and Jarić-Ebenhard do not use the observation origin at all.

Endpoint estimators

robson1964() uses the final gap for a Quenouille-type bias-reduced point estimate and a one-sided confidence interval. strauss1989() supplies an unbiased endpoint estimate and one-sided interval under uniform occurrence. ole() fits the Weibull extreme-value model to the k most recent sighting events and reports the two-sided interval of Roberts & Solow (2003).

robson1964(sd)
#> <Robson & Whitlock (1964)>
#>   estimate: 1939
#>   95% one-sided CI: [1936, 1993]
strauss1989(sd)
#> <Strauss & Sadler (1989)>
#>   estimate: 1938.25
#>   95% one-sided CI: [1936, 1943.41]
ole(sd)
#> <OLE (Roberts & Solow 2003)>
#>   estimate: 1941.97
#>   95% CI: [1937.48, 1956.75]
ole(sd, k = 10)
#> <OLE (Roberts & Solow 2003)>
#>   estimate: 1940.43
#>   95% CI: [1936.35, 1957.28]

Persistence tests

solow1993(sd, test_year = 2000)
#> <Solow (1993a)>
#>   95% upper confidence bound: 1939.32
solow1993b(sd, test_year = 2000)
#> Warning: p-value never falls to alpha before `test_year`; returning NA.
#> <Solow (1993b)>
#>   estimate: 1940.97
#>   95% upper confidence bound: NA
solow2005(sd, test_year = 2000)
#> <Solow (2005) Weibull test>
#>   95% upper confidence bound: 1953.99
mcinerny2006(sd, test_year = 2000)
#> <McInerny, Roberts, Davy & Cribb (2006)>
#>   95% upper confidence bound: 1942
burgman1995(sd, test_year = 2000)
#> <Burgman, Grimson & Ferson (1995), equation 2>
#>   95% upper confidence bound: 1940
solow_roberts2003(sd, test_year = 2000)
#> <Solow & Roberts (2003)>
#>   95% upper confidence bound: 1993
jaric2010(sd, test_year = 2000)
#> <Jaric & Ebenhard (2010)>
#>   95% upper confidence bound: 1980.02

solow1993() assumes a stationary Poisson process.
solow1993b() allows an exponentially declining sighting rate. Its upper bound may not exist: the p-value has a positive floor however late the candidate time.
solow2005() uses the Weibull extreme-value p-value associated with OLE. Its bound is one-sided, so it equals the upper limit of ole() at twice the significance level.
mcinerny2006() conditions on the previous binary sighting rate.
burgman1995() uses discrete periods. By default it applies equation 2, the test for extinction based on the final run of empty periods; equation = 4 tests for a change in conservation status using the longest empty run anywhere in the record.
solow_roberts2003() uses only the two most recent distinct sightings; its bound coincides with that of robson1964().
jaric2010() uses the average interval and can adjust it for a trend in consecutive interval lengths. It is an index rather than a formal test.

The full p-value curve is available with data_out = TRUE:

curve <- jaric2010(sd, test_year = 2000, data_out = TRUE)
plot(curve$time, curve$chance, type = "l",
     xlab = "candidate time", ylab = "index")
abline(h = 0.05, lty = 2)

A Bayesian view

Under the constant-rate model, solow1993_bayes() gives the Bayes factor of Solow (1993a) in favour of persistence to a candidate time, assuming a uniform prior on the extinction time and a noninformative prior on the sighting rate. Combined with a prior probability of persistence, it gives the posterior probability that the species is extant, which is the quantity the p-values above do not provide.

solow1993_bayes(sd, test_year = 2000, prior = 0.5)
#> <Solow (1993a) Bayes factor>
#>   Bayes factor for persistence to 2000: 7.525e-14
#>   posterior probability of persistence (prior 0.5): 7.525e-14

Uncertain sightings

jaric_roberts2014() assigns a probability of validity to every sighting. EDE requires binary counts because reliability belongs to individual records; this is a package restriction, not one of the source paper. A sighting at the observation origin fixes the origin and its reliability is not used.

uncertain <- sighting_data(data.frame(
  year = c(1900, 1910, 1920, 1930),
  sightings = 1
))

jaric_roberts2014(
  uncertain,
  reliability = c(1.0, 0.9, 0.6, 0.3)
)
#> <Jaric & Roberts (2014)>
#>   estimate: 1930.99
#>   95% upper confidence bound: 2005.22

The result includes the effective number of sightings and the reliability-adjusted endpoint as additional components. Supply test_year to add the reliability-adjusted p-value of equations 7-8. Reliabilities should be elicited independently of the extinction analysis.

Other approaches to uncertain sightings:

thompson2013() handles several observation types (for example specimens, expert sightings and local reports), each with its own range of detection and validity probabilities, in discrete time.
brook2019() resamples the record, keeping each sighting with probability equal to its reliability, and applies any EDE estimator to every sub-record. Solow (2023) criticised the statistical properties of the resulting date; the authors present it as pre-processing rather than an estimator, and its help page discusses both views.

Other Bayesian models

rout2009() generalises the Solow (1993a) Bayes factor to a sighting rate that declines as a power of time, and rout2009_absences() turns it into a decision rule: the number of consecutive absent surveys after which declaring eradication minimises the expected cost. caley_barry2014() models yearly detections with constant detection and extinction probabilities; its posterior is computed exactly:

fox <- sighting_data(data.frame(year = c(2001, 2003, 2005, 2006), count = 1))
caley_barry2014(fox, test_year = 2012)
#> <Caley & Barry (2014) constant model (published)>
#>   Bayes factor for persistence to 2012: 0.9147
#>   posterior probability of persistence (prior 0.0680721): 0.06263

Records and surveys

thompson2017() implements the iterative model of Thompson et al. (2017), which underlies the IUCN Red List tool for estimating the probability that a taxon is extinct. Each year contributes either a record, with a probability of being correct, or a survey that failed to find the species, described by its extent, its chance of detecting the species and its reliability. Ranges for these quantities give interval and Monte Carlo bounds. thompson2019() casts the same data as Bayesian updating, in which each year’s posterior is the next year’s prior.

Fossil ranges

The same question arises for the last appearance of a taxon in the fossil record. strauss1989() assumes fossil horizons are uniformly distributed; marshall1994() drops that assumption and bounds the confidence interval using only the distribution of gap sizes, with marshall1994_trend() to test whether gap size changes through the section; marshall1997() accepts a fossil recovery potential that varies through the section; and weiss_marshall1999() gives a Bayesian estimate for sections sampled at discrete horizons, with absences recorded. EDE always estimates the end-point at the largest time, so ages in millions of years should be entered as time = -age, and first appearances as negative stratigraphic positions.

Comparing methods

ede_compare() applies every method that needs only the sighting record:

ede_compare(sighting_data(monk_seal), test_year = 2050)
#>                                          method     function_name estimate
#> 1                    OLE (Roberts & Solow 2003)               ole 1963.528
#> 2                      Robson & Whitlock (1964)        robson1964 1956.000
#> 3                       Strauss & Sadler (1989)       strauss1989 1961.250
#> 4                                 Solow (1993a)         solow1993       NA
#> 5                                 Solow (1993b)        solow1993b 1957.605
#> 6                     Solow (2005) Weibull test         solow2005       NA
#> 7        McInerny, Roberts, Davy & Cribb (2006)      mcinerny2006       NA
#> 8  Burgman, Grimson & Ferson (1995), equation 2       burgman1995       NA
#> 9                        Solow & Roberts (2003) solow_roberts2003       NA
#> 10                      Jaric & Ebenhard (2010)         jaric2010       NA
#>      lower    upper interval_type
#> 1  1952.67 2093.799     two-sided
#> 2  1952.00 2028.000     one-sided
#> 3  1952.00 1993.245     one-sided
#> 4       NA 1993.245     one-sided
#> 5       NA       NA     one-sided
#> 6       NA 2029.444     one-sided
#> 7       NA 1979.000     one-sided
#> 8       NA 1994.000     one-sided
#> 9       NA 2028.000     one-sided
#> 10      NA       NA     one-sided
#>                                                              note
#> 1                                                                
#> 2                                                                
#> 3                                                                
#> 4                                                                
#> 5  p-value never falls to alpha before `test_year`; returning NA.
#> 6                                                                
#> 7                                                                
#> 8                                                                
#> 9                                                                
#> 10 p-value never falls to alpha before `test_year`; returning NA.

Choosing a method

Use OLE when the recent record is sufficiently long and an extreme-value model is defensible.
Use solow1993() only when a constant sighting rate is plausible.
Use solow1993b() for an approximately exponential decline.
Use jaric2010() when the intervals themselves show a gradual trend.
Use burgman1995() for discrete frequency data; use equation = 4 only for questions about a change in status.
Use jaric_roberts2014() when observations have explicit reliability assessments, and thompson2013() when they fall into types with uncertain detection probabilities.
Use thompson2017() when unsuccessful surveys are documented, as in an IUCN Red List assessment.
Use the Bayesian models when a posterior probability that the species is extant is needed; report the prior alongside it.

Running several scientifically defensible methods is useful, but p-values from models whose assumptions are violated should not be combined or interpreted as extinction probabilities.

References

References for the methods introduced after the core estimators (Rout et al. 2009; Thompson et al. 2013, 2017, 2019; Caley & Barry 2014; Brook et al. 2019; Marshall 1994, 1997; Weiss & Marshall 1999) are given on the help page of each function.

Burgman, M. A., Grimson, R. C., & Ferson, S. (1995). Inferring threat from scientific collections. Conservation Biology, 9(4), 923-928.

Cooke, P. (1980). Optimal linear estimation of bounds of random variables. Biometrika, 67, 257-258.

Jarić, I., & Ebenhard, T. (2010). A method for inferring extinction based on sighting records that change in frequency over time. Wildlife Biology, 16(3), 267-275.

Jarić, I., & Roberts, D. L. (2014). Accounting for observation reliability when inferring extinction based on sighting records. Biodiversity and Conservation, 23(11), 2801-2815.

McInerny, G. J., Roberts, D. L., Davy, A. J., & Cribb, P. J. (2006). Significance of sighting rate in inferring extinction and threat. Conservation Biology, 20(2), 562-567.

Roberts, D. L., & Solow, A. R. (2003). When did the dodo become extinct? Nature, 426, 245.

Robson, D. S., & Whitlock, J. H. (1964). Estimation of a truncation point. Biometrika, 51, 33-39.

Solow, A. R. (1993a). Inferring extinction from sighting data. Ecology, 74, 962-964.

Solow, A. R. (1993b). Inferring extinction in a declining population. Journal of Mathematical Biology, 32, 79-82.

Solow, A. R. (2005). Inferring extinction from a sighting record. Mathematical Biosciences, 195, 47-55.

Solow, A. R., & Roberts, D. L. (2003). A nonparametric test for extinction based on a sighting record. Ecology, 84, 1329-1332.

Strauss, D., & Sadler, P. M. (1989). Classical confidence intervals and Bayesian probability estimates for ends of local taxon ranges. Mathematical Geology, 21, 411-427.