Title:	Extinction Risk Estimation
Version:	1.1.0
Description:	Estimates extinction risk from population time series under a drifted Wiener process using MLE and observation-error-and-autocovariance-robust estimators, with confidence intervals based on the w-z method.
License:	BSD_2_clause + file LICENSE
Encoding:	UTF-8
Language:	en-US
RoxygenNote:	7.3.3
NeedsCompilation:	no
Packaged:	2026-03-17 07:58:09 UTC; hako
Author:	Hiroshi Hakoyama [aut, cre, cph]
Maintainer:	Hiroshi Hakoyama <hiroshi.hakoyama@gmail.com>
Repository:	CRAN
Date/Publication:	2026-03-17 08:10:13 UTC

extr: Extinction Risk Estimation

Description

Estimates extinction risk from population time series under a drifted Wiener process using MLE and observation-error-and-autocovariance-robust estimators, with confidence intervals based on the w-z method.

Author(s)

Maintainer: Hiroshi Hakoyama hiroshi.hakoyama@gmail.com (ORCID) [copyright holder]

Extinction Risk Estimation for a Density-Independent Model

Description

Estimates demographic parameters and extinction probability under a density-independent (drifted Wiener) model. From a time series of population sizes, it estimates growth rate and variance using either a naive MLE variance estimator or an observation-error-and-autocovariance-robust (OEAR) estimator, then evaluates extinction risk over a horizon t^{\ast} with confidence intervals from the w–z method. In the observation-error-free setting, the naive w–z interval achieves near-nominal coverage, while OEAR can be more robust under additive observation error and unknown short-run autocovariance.

Usage

ext_di(
  dat,
  ne = 1,
  th = 100,
  alpha = 0.05,
  unit = "years",
  qq_plot = FALSE,
  formatted = TRUE,
  digits = getOption("extr.digits", 5L),
  method = c("naive", "oear")
)

Arguments

Details

Population dynamics follow

dX=\mu\,dt+\sigma\,dW,

where X(t)=\log N(t) and W is a Wiener process. Extinction risk is

G=\Pr[T\le t^{\ast}\mid N(0)=n_0,n_e,\mu,\sigma],

the probability that population size falls below n_e within t^{\ast}. Irregular observation intervals are allowed.

The function first estimates \mu, then estimates a variance quantity according to method:

1. naive: uses the drifted-Wiener MLE variance (Dennis et al., 1991).

2. oear: uses an observation-error-and-autocovariance-robust (OEAR) HAC long-run-variance estimator with AR(1) pre-whitening and a Bartlett kernel (Hakoyama, in press).

For method = "oear", robustness to additive observation error relies on long-run-variance cancellation of differenced observation error (McNamara and Harding, 2004), while robustness to unknown short-run autocovariance comes from HAC estimation of long-run variance.

Extinction probability is then computed as G(w,z) (Lande and Orzack, 1988), and confidence intervals are constructed using the w–z method (Hakoyama, in press).

Value

A list (class "ext_di" when formatted=TRUE) containing extinction-risk estimates and diagnostics. Core elements include Growth.rate, Variance, Unbiased.variance, lower_cl_s, upper_cl_s, linear_g, log_g, log_q, ci_linear_g, ci_log_g, ci_log_q, and data/input summaries (nq, xd, sample.size, unit, ne, th, alpha).

Variance is method-dependent (naive MLE for method="naive", OEAR HAC estimate for method="oear"). aic is reported only for method="naive" and is NA for method="oear". method_diag is NULL for naive and contains rho_tilde_pw and j for oear. lower_cl_s and upper_cl_s are chi-square based: exact under method="naive" model assumptions, and pragmatic plug-in approximations for method="oear".

References

Andrews, D. W. K. (1991). Heteroskedasticity and autocorrelation consistent covariance matrix estimation. Econometrica, 59(3), 817-858.

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. Confidence intervals for extinction risk: validating population viability analysis with limited data. Methods in Ecology and Evolution. In press. Preprint, doi:10.48550/arXiv.2509.09965

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.

McNamara, J. M. and Harding, K. C. (2004). Measurement error and estimates of population extinction risk. Ecology Letters, 7(1), 16-20.

Newey, W. K. and West, K. D. (1987). A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica, 55(3), 703-708.

See Also

statistics_di, oear_sigma2_hac, extinction_probability_di, confidence_interval_wz_di, print.ext_di

Examples

# Example from Dennis et al. (1991), Yellowstone grizzly bears.
# Population is a running 3-year sum (3-year moving total) digitized from
# Fig. 5.
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)

# Probability of decline to 10 individuals within 100 years
ext_di(dat, th = 100, ne = 10)

# With QQ-plot
ext_di(dat, th = 100, ne = 10, qq_plot = TRUE)

# OEAR method
ext_di(dat, th = 100, method = "oear")

Computes Extinction Probability for a Density-Independent Model

Description

Compute G(w,z) and its complement Q(w,z)=1-G(w,z) for a density-independent (drifted Wiener) model, on both linear and log scales.

Usage

ext_prob_di(w, z)

log_ext_prob_di(w, z)

log_ext_comp_di(w, z)

ext_prob_format_di(w, z, digits = 5L)

Arguments

Details

For any t>0 with w+z>0,

\Pr[T \leq t] = G(w,z)=\Phi(-w)+ \exp\!\left(\tfrac{z^2-w^2}{2}\right)\Phi(-z), \qquad Q(w,z)=1-G(w,z).

Here \Phi and \phi denote the standard normal CDF and PDF.

Stability strategy. (i) For large z, rewrite the product \exp((z^2-w^2)/2)\,\Phi(-z) via the Mills ratio and replace it by an 8-term asymptotic series when z \ge 19. (ii) On the log scale, use log-sum-exp and a stable log-difference (log1mexp) built from log1p/expm1 to retain tail info.

Domain. Scalar inputs are assumed and require w+z>0.

Functions.

• ext_prob_di(w,z): returns G(w,z) (linear scale).

• log_ext_prob_di(w,z): returns \log G(w,z).

• log_ext_comp_di(w,z): returns \log Q(w,z).

• ext_prob_format_di(w,z,digits): formats a point estimate using repr_mode() and format_by_mode().

Value

For ext_prob_di: numeric scalar G(w,z).
For log_ext_prob_di: numeric scalar \log G(w,z).
For log_ext_comp_di: numeric scalar \log Q(w,z).
For ext_prob_format_di: character string formatted for display.

Author(s)

Hiroshi Hakoyama, hiroshi.hakoyama@gmail.com

See Also

statistics_di, repr_mode, format_by_mode

Format a Probability using a Display Mode

Description

Format a probability (point estimate or CI bound) according to a chosen representation mode. Allowed modes are "linear_G", "log_G", "log_Q", and "undefined". The function performs no inference; it only formats the numeric values already computed (linear and log).

Notes on arguments:

• log_g, log_q, ci_log_g, ci_log_q are natural logarithms (from log() in R).

• For "log_Q": the upper bound for G uses ci_log_q[1] and the lower bound uses ci_log_q[2]. This ordering reflects that Q = 1-G decreases as G increases.

Usage

format_by_mode(
  mode,
  kind,
  linear_g,
  log_g,
  log_q,
  ci_linear_g,
  ci_log_g,
  ci_log_q,
  digits
)

Arguments

Value

A character scalar formatted for display. Returns "NA" if the required value is non-finite or missing.

Author(s)

Hiroshi Hakoyama, hiroshi.hakoyama@gmail.com

OEAR Diffusion-Scale Estimator

Description

Compute the OEAR (observation-error-and-autocovariance-robust) diffusion-scale estimate \tilde{\sigma}^2 by estimating the long-run variance (LRV) of standardized growth increments using AR(1) pre-whitening and a Bartlett (Newey-West) HAC estimator, followed by recoloring.

Arguments

Details

The estimator targets the LRV of the standardized centered increment series

U_i = (\Delta Y_i - \hat{\mu}\,\tau_i)/\sqrt{\tau_i}, \qquad i=1,\ldots,q,

which is interpreted as an effective environmental variance per unit time.

Implementation steps:

1. Center u_i by subtracting \bar u to obtain \tilde u_i.

2. Estimate the AR(1) pre-whitening coefficient \tilde\rho_{\mathrm{pw}} by OLS in \tilde u_i = \rho_{\mathrm{pw}}\,\tilde u_{i-1} + \varepsilon_i and form pre-whitened residuals \tilde\varepsilon_i.

3. Compute residual autocovariances and the Bartlett HAC LRV \widetilde{\mathcal C}^{(\varepsilon)}_{\mathrm{NW}}(J).

4. Choose the truncation lag J via the Andrews (1991) AR(1) plug-in rule specialized to the Bartlett window.

5. Recolor to the original scale to obtain

   \tilde{\sigma}^2 = \widetilde{\mathcal C}_{\mathrm{NW}}(J) = \widetilde{\mathcal C}^{(\varepsilon)}_{\mathrm{NW}}(J)/ (1-\tilde\rho_{\mathrm{pw}})^2.

This construction is designed to be robust to short-run autocovariance and, under additive observation error on the log scale, to the cancellation property of the LRV when increments are appropriately centered.

Value

A list with components:

sigma2_tilde

Numeric scalar. OEAR diffusion-scale estimate \tilde{\sigma}^2.

rho_tilde_pw

Numeric scalar. Estimated AR(1) pre-whitening coefficient \tilde\rho_{\mathrm{pw}}.

j

Integer. Selected Bartlett truncation lag $J$.

References

Andrews, D. W. K. (1991). Heteroskedasticity and autocorrelation consistent covariance matrix estimation. Econometrica, 59(3), 817–858.

Newey, W. K. and West, K. D. (1987). A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica, 55(3), 703–708.

Print Extinction Risk Estimates

Description

S3 method that prints a formatted summary for an ext_di object. The output includes the plug-in extinction probability and its CI, method-specific parameter estimates, a brief data summary, and input settings.

Usage

## S3 method for class 'ext_di'
print(x, digits = NULL, ...)

Arguments

Value

Invisibly returns x after printing the formatted results.

Display Modes for Extinction Probabilities

Description

A pair of internal functions for formatting extinction probabilities: repr_mode() decides how to represent a probability G based on its magnitude, and format_by_mode() formats numbers accordingly for presentation. These functions affect display only and do not alter underlying calculations.

Usage

repr_mode(g)

Arguments

Details

Probabilities very close to 0 or 1 are displayed on a log scale to improve readability. The thresholds are fixed at 1e-300 ("near zero") and 1 - 1e-15 ("near one"). Non-finite inputs are reported as "undefined".

Value

repr_mode() returns a character scalar, one of:

• "linear_G" – display G on the linear scale;

• "log_G" – display G on log10 scale for tiny probabilities;

• "log_Q" – display Q=1-G on log10 scale when G \approx 1;

• "undefined" – input was not finite.

format_by_mode() returns a character scalar formatted according to the selected mode and the kind of value (point estimate, lower, or upper CI bound).

Author(s)

Hiroshi Hakoyama, hiroshi.hakoyama@gmail.com

Computes the w and z Statistics

Description

Estimators for the w and z statistics used in extinction probability calculations under a density-independent population model.

Usage

w_statistic(mu, xd, s, th)

z_statistic(mu, xd, s, th)

Arguments

Details

The statistics are defined as

\hat w = \frac{\hat \mu t^{\ast} + x_d}{\sqrt{\hat \sigma^2 t^{\ast}}}, \qquad \hat z = \frac{- \hat \mu t^{\ast} + x_d}{\sqrt{\hat \sigma^2 t^{\ast}}}.

Value

numeric: Value of the statistic.

Author(s)

Hiroshi Hakoyama, hiroshi.hakoyama@gmail.com

Confidence Intervals for Extinction Probability (w–z, DI model)

Description

Computes confidence intervals (CIs) for extinction probability in a density-independent (drifted Wiener) model using the w–z method, and provides a formatter for display-ready CI strings.

Usage

confidence_interval_wz_di(mu, xd, s, th, tq, qq, alpha, prob_fun = ext_prob_di)

ci_wz_format_di(mu, xd, s, th, tq, qq, alpha, digits = 5L)

Arguments

Details

The w–z method derives CIs by inverting noncentral-t distributions for the transformed statistics w and z, then combining them to form bounds on G(w,z).

Exact confidence intervals for w and z are obtained by numerically solving the noncentral-t quantile equations corresponding to the observed statistics.

The CI for G(w,z) is then approximated by

\bigl( G(\overline{w},\,\underline{z}),\; G(\underline{w},\,\overline{z}) \bigr),

where \overline{w}, \underline{w}, \overline{z}, and \underline{z} are the upper and lower confidence limits for w and z, respectively. Across the full parameter space, this approach achieves near-nominal coverage. When $z$ is large and positive, $G$ depends only on $w$, so exact CIs are available.

The argument prob_fun selects which probability is evaluated:

• ext_prob_di, log_ext_prob_di: CIs for G(w,z) and \log G(w,z); returned as (lower, upper).

• log_ext_comp_di: CIs for \log Q(w,z)= \log (1-G(w,z)); returned as (lower, upper).

Value

For confidence_interval_wz_di: numeric vector c(lower, upper) on the chosen scale (natural log if log_*).
For ci_wz_format_di: named character vector c(lower, upper) with values preformatted for display.

Author(s)

Hiroshi Hakoyama, hiroshi.hakoyama@gmail.com