Extinction Risk Estimation for a Density-Independent Model
ext_di.Rd
Estimates demographic parameters and extinction probability under a density-independent (drifted Wiener) model. From a time series of population sizes, it computes MLEs of growth rate and environmental variance, then evaluates extinction risk over a horizon \(t^{\ast}\). Confidence intervals are constructed by the \(w\)-\(z\) method, which achieve near-nominal coverage across the full parameter space.
Usage
ext_di(
dat,
ne = 1,
th = 100,
alpha = 0.05,
unit = "years",
qq_plot = FALSE,
formatted = TRUE,
digits = getOption("extr.digits", 5L)
)
Arguments
- dat
Data frame with two numeric columns: time (strictly increasing) and population size. Column names are not restricted.
- ne
Numeric. Extinction threshold \(n_e \ge 1\). Default is 1.
- th
Numeric. Time horizon \(t^{\ast} > 0\). Default is 100.
- alpha
Numeric. Significance level \(\alpha \in (0,1)\). Default is 0.05.
- unit
Character. Unit of time (e.g., "years", "days", "generations"). Default is "years".
- qq_plot
Logical. If
TRUE
, draws a QQ-plot of standardized increments to check model assumptions. Default isFALSE
.- formatted
Logical. If
TRUE
, returns an"ext_di"
object; otherwise returns a raw list. Default isTRUE
.- digits
Integer. Preferred significant digits for printing. Affects display only. Default is
getOption("extr.digits", 5)
.
Value
A list (class "ext_di"
if formatted=TRUE
) with:
Growth.rate
,Variance
,Unbiased.variance
;AIC
;Extinction.probability
with confidence limits;data summary (
nq
,xd
,sample.size
);input parameters (
unit
,ne
,th
,alpha
).
Details
Population dynamics follow $$dX=\mu\,dt+\sigma\,dW,$$ where \(X(t)=\log N(t)\), \(\mu\) is the growth rate, \(\sigma^2\) the environmental variance, and \(W\) a Wiener process. Extinction risk is $$G=\Pr[T\le t^{\ast}\mid N(0)=n_0,n_e,\mu,\sigma],$$ the probability the population falls below \(n_e\) within \(t^{\ast}\). Irregular intervals are allowed.
The function:
estimates \(\mu\) and \(\sigma^2\) (Dennis et al., 1991),
computes extinction probability \(G(w,z)\) (Lande and Orzack, 1988),
constructs confidence intervals for \(G\) using the \(w\)-\(z\) method (Hakoyama, 2025).
Numerical range. Probabilities are evaluated on \(G\), \(\log G\), and
\(\log(1-G)\) scales. The log-scale removes the \(\approx
4.94\times10^{-324}\) lower bound of linear doubles and extends the safe
range down to exp(-DBL_MAX)
(kept symbolically), avoiding
underflow/cancellation.
References
Lande, R. and Orzack, S.H. (1988) Extinction dynamics of age-structured populations in a fluctuating environment. Proceedings of the National Academy of Sciences, 85(19), 7418–7421.
Dennis, B., Munholland, P.L., and Scott, J.M. (1991) Estimation of growth and extinction parameters for endangered species. Ecological Monographs, 61, 115–143.
Hakoyama, H. (2025) Confidence intervals for extinction risk: validating population viability analysis with limited data. Preprint, doi:10.48550/arXiv.2509.09965
Author
Hiroshi Hakoyama, hiroshi.hakoyama@gmail.com
Examples
# Example from Dennis et al. (1991), Yellowstone grizzly bears
dat <- data.frame(Time = 1959:1987,
Population = c(44, 47, 46, 44, 46, 45, 46, 40, 39, 39, 42, 44, 41, 40,
33, 36, 34, 39, 35, 34, 38, 36, 37, 41, 39, 51, 47, 57, 47))
# Probability of decline to 1 individual within 100 years
ext_di(dat, th = 100)
#> --- Estimates ---
#> Estimate
#> Probability of decline to 1 within 100 years (MLE): 9.4128e-05
#> Growth rate (MLE): 0.0023556
#> Variance (MLE): 0.01087
#> Unbiased variance: 0.011273
#> AIC for the distribution of N: 165.06
#> CI
#> Probability of decline to 1 within 100 years (MLE): (1.4586e-13, 0.5653)
#> Growth rate (MLE): (-0.038814, 0.043525)
#> Variance (MLE): (0.0070464, 0.020885)
#> Unbiased variance: -
#> AIC for the distribution of N: -
#>
#> --- Data Summary ---
#> Value
#> Current population size, nq: 47
#> xd = ln(nq / ne): 3.8501
#> Sample size, q + 1: 29
#>
#> --- Input Parameters ---
#> Parameter
#> Time unit: years
#> Extinction threshold of population size, ne: 1
#> Time window for extinction risk evaluation (years), th: 100.0
#> Significance level, alpha: 0.05
# Probability of decline to 10 individuals within 100 years
ext_di(dat, th = 100, ne = 10)
#> --- Estimates ---
#> Estimate
#> Probability of decline to 10 within 100 years (MLE): 0.096852
#> Growth rate (MLE): 0.0023556
#> Variance (MLE): 0.01087
#> Unbiased variance: 0.011273
#> AIC for the distribution of N: 165.06
#> CI
#> Probability of decline to 10 within 100 years (MLE): (1.0699e-05, 0.9898)
#> Growth rate (MLE): (-0.038814, 0.043525)
#> Variance (MLE): (0.0070464, 0.020885)
#> Unbiased variance: -
#> AIC for the distribution of N: -
#>
#> --- Data Summary ---
#> Value
#> Current population size, nq: 47
#> xd = ln(nq / ne): 1.5476
#> Sample size, q + 1: 29
#>
#> --- Input Parameters ---
#> Parameter
#> Time unit: years
#> Extinction threshold of population size, ne: 10
#> Time window for extinction risk evaluation (years), th: 100.0
#> Significance level, alpha: 0.05
# With QQ-plot
ext_di(dat, th = 100, qq_plot = TRUE)
#> --- Estimates ---
#> Estimate
#> Probability of decline to 1 within 100 years (MLE): 9.4128e-05
#> Growth rate (MLE): 0.0023556
#> Variance (MLE): 0.01087
#> Unbiased variance: 0.011273
#> AIC for the distribution of N: 165.06
#> CI
#> Probability of decline to 1 within 100 years (MLE): (1.4586e-13, 0.5653)
#> Growth rate (MLE): (-0.038814, 0.043525)
#> Variance (MLE): (0.0070464, 0.020885)
#> Unbiased variance: -
#> AIC for the distribution of N: -
#>
#> --- Data Summary ---
#> Value
#> Current population size, nq: 47
#> xd = ln(nq / ne): 3.8501
#> Sample size, q + 1: 29
#>
#> --- Input Parameters ---
#> Parameter
#> Time unit: years
#> Extinction threshold of population size, ne: 1
#> Time window for extinction risk evaluation (years), th: 100.0
#> Significance level, alpha: 0.05
# Change digits
ext_di(dat, th = 100, ne = 10, digits = 9)
#> --- Estimates ---
#> Estimate
#> Probability of decline to 10 within 100 years (MLE): 0.0968521443
#> Growth rate (MLE): 0.00235564171
#> Variance (MLE): 0.0108702457
#> Unbiased variance: 0.0112728474
#> AIC for the distribution of N: 165.058678
#> CI
#> Probability of decline to 10 within 100 years (MLE): (1.0698563e-05, 0.989801572)
#> Growth rate (MLE): (-0.0388142081, 0.0435254915)
#> Variance (MLE): (0.00704642493, 0.0208851222)
#> Unbiased variance: -
#> AIC for the distribution of N: -
#>
#> --- Data Summary ---
#> Value
#> Current population size, nq: 47
#> xd = ln(nq / ne): 1.5476
#> Sample size, q + 1: 29
#>
#> --- Input Parameters ---
#> Parameter
#> Time unit: years
#> Extinction threshold of population size, ne: 10
#> Time window for extinction risk evaluation (years), th: 100.0
#> Significance level, alpha: 0.05