Mark-Recapture and Exploitation Rate Estimation

Chapman estimator (default)

The Chapman estimator is a bias-corrected Petersen estimator recommended when the recapture count is small. With \(M\) tagged anglers released, \(n\) anglers checked in the second sample, and \(m\) recaptures:

\[\hat{N} = \frac{(M+1)(n+1)}{(m+1)} - 1\]

Variance:

\[\widehat{\text{Var}}(\hat{N}) = \frac{(M+1)(n+1)(M-m)(n-m)}{(m+1)^2(m+2)}\]

library(tidycreel)

result_chapman <- estimate_angler_n(M = 200L, n = 50L, m = 10L)
print(result_chapman)
#> 
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: mark-recapture-chapman
#> Variance: chapman
#> Confidence level: 95%
#> 
#> # A tibble: 1 × 6
#>   parameter estimate    se ci_lower ci_upper     n
#>   <chr>        <dbl> <dbl>    <dbl>    <dbl> <int>
#> 1 N_hat         931.  232.     605.    1715.    10

The n column in the estimates tibble records the recapture count (m), which determines precision.

Petersen estimator

The unadjusted Lincoln–Petersen estimator:

\[\hat{N} = \frac{M \cdot n}{m}\]

tidycreel enforces a minimum of \(m \geq 7\) recaptures; below this threshold the Petersen estimator carries large positive bias and Chapman should be used instead.

result_petersen <- estimate_angler_n(M = 200L, n = 50L, m = 10L, method = "petersen")
print(result_petersen)
#> 
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: mark-recapture-petersen
#> Variance: petersen
#> Confidence level: 95%
#> 
#> # A tibble: 1 × 6
#>   parameter estimate    se ci_lower ci_upper     n
#>   <chr>        <dbl> <dbl>    <dbl>    <dbl> <int>
#> 1 N_hat         1000  283.     605.    1715.    10

Schnabel estimator (multi-occasion)

For \(K \geq 2\) successive sampling occasions, the Schnabel estimator pools information across occasions, carrying Chapman’s (1952) small-sample correction to the recapture total:

\[\hat{N} = \frac{\sum_{k=1}^{K} M_k n_k}{\sum_{k=1}^{K} m_k + 1}\]

where \(M_k\) is the cumulative count of marked-at-large anglers before occasion \(k\) (so \(M_1 = 0\)), \(n_k\) is the catch on occasion \(k\), and \(m_k\) is the number of recaptures.

The \(+1\) is the default as of tidycreel 3.0.0. Each \(m_k\) is approximately Poisson with parameter \(M_k n_k / N\), which is what motivates the correction. Dettloff (2023) simulated both forms and found the unadjusted estimator turns biased high at moderate sample sizes — inflating harvest downstream — while the adjusted form’s bias approaches zero without ever becoming positive, at lower variance and no cost in large samples. It also keeps bias handling consistent across occasion counts: Schnabel reduces to Lincoln-Petersen at \(K = 2\), so an unadjusted Schnabel would have silently dropped the correction that method = "chapman" applies to the same data. Pass bias_adjust = FALSE for the pre-3.0.0 form, which is also what fishmethods::schnabel() computes.

The correction is worth \(-1/(\sum m_k + 1)\) in relative terms: \(-33\%\) at \(\sum m_k = 2\), \(-1.9\%\) at 52, and \(-0.2\%\) at 500. It matters at low recapture totals and is negligible at high ones.

Confidence intervals use the Poisson distribution when \(\sum m_k < 50\) and a normal approximation on \(1/\hat{N}\) otherwise. The Poisson branch inverts the distribution of \(\sum m_k\) directly rather than centring on \(\hat{N}\), so it is a valid interval under either form and its bounds do not move with bias_adjust; the normal branch is built around \(1/\hat{N}\) and does.

result_schnabel <- estimate_angler_n(
  M      = c(0L, 47L, 91L, 131L),
  n      = c(50L, 50L, 50L, 50L),
  m      = c(0L,  4L,  6L,  8L),
  method = "schnabel"
)
print(result_schnabel)
#> 
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: mark-recapture-schnabel
#> Variance: delta
#> Confidence level: 95%
#> 
#> # A tibble: 1 × 6
#>   parameter estimate    se ci_lower ci_upper     n
#>   <chr>        <dbl> <dbl>    <dbl>    <dbl> <int>
#> 1 N_hat         708.  158.     498.     1345    18

Four sampling occasions yield 18 total recaptures here, so the Poisson CI branch is used automatically.

Schumacher-Eschmeyer estimator (multi-occasion regression)

The same data can be pooled by regression instead of by ratio. Schumacher and Eschmeyer (1943) fit \(m_k/n_k\) against \(M_k\) through the origin with slope \(1/N\), which Seber (1982, sec. 4.1.3) develops as

\[\hat{N} = \frac{\sum_{k} n_k M_k^2}{\sum_{k} m_k M_k}\]

with the interval from his eq. (4.17) on \(K - 2\) degrees of freedom. The degrees of freedom differ from Schnabel’s \(K - 1\) for a stated reason: \(y_1\) is identically zero when \(M_1 = 0\), so it is not a random observation and does not enter the regression. It needs at least three occasions.

result_se <- estimate_angler_n(
  M      = c(0L, 47L, 91L, 131L),
  n      = c(50L, 50L, 50L, 50L),
  m      = c(0L,  4L,  6L,  8L),
  method = "schumacher"
)
print(result_se)
#> 
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: mark-recapture-schumacher
#> Variance: delta
#> Confidence level: 95%
#> 
#> # A tibble: 1 × 6
#>   parameter estimate    se ci_lower ci_upper     n
#>   <chr>        <dbl> <dbl>    <dbl>    <dbl> <int>
#> 1 N_hat         699.  44.7     548.     964.    18

Which of the two to use is a design question. Seber expects the regression form “to be robust with regard to departures from the underlying assumptions” and recommends it “in conjunction with the other methods” — as a cross-check rather than a replacement. Dettloff (2023) found the two adjusted forms effectively equivalent at larger sample sizes, with Schumacher-Eschmeyer less variable and Schnabel reaching unbiasedness slightly sooner.

Do not choose whichever gives the narrower interval. Hansen and Van Kirk (2018) computed both and selected the estimator producing the smallest 95% CI, and that procedure does not have 95% coverage — picking the narrower of two intervals after seeing them conditions on the luckier draw. tidycreel does not implement the rule. Decide before looking, or report both.

Unstratified

\[\hat{u} = \frac{C \cdot m}{T \cdot n}, \quad \widehat{\text{Var}}(\hat{u}) \approx \left(\frac{C}{T}\right)^2 \frac{p(1-p)}{n} + p^2 \frac{SE_C^2}{T^2}\]

where \(p = m/n\).

result_expl <- estimate_exploitation_rate(
  T    = 200L,
  C    = 450.0,
  se_C = 42.0,
  n    = 180L,
  m    = 15L
)
#> ℹ `C` was supplied as a bare number, so its target cannot be verified.
#>   It must be a season (period) total -- a sampled-day total understates `u` by
#>   the sampling fraction.
#>   Pass the `estimate_total_harvest()` result itself to have this checked.
print(result_expl)
#> 
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: Exploitation Rate (Mark-Recapture)
#> Variance: delta
#> Confidence level: 95%
#> Unit: proportion
#> 
#> # A tibble: 1 × 8
#>   estimate     se ci_lower ci_upper     n     T     C     m
#>      <dbl>  <dbl>    <dbl>    <dbl> <int> <int> <dbl> <int>
#> 1    0.188 0.0495   0.0897    0.285   180   200   450    15

Stratified (T-weighted)

When the survey spans multiple strata (day types, access areas), supply a data frame with one row per stratum. The T-weighted aggregate is automatically appended as an .overall row:

\[\hat{u}_{overall} = \frac{\sum_h T_h \hat{u}_h}{\sum_h T_h}, \quad \widehat{\text{Var}}(\hat{u}_{overall}) = \frac{\sum_h T_h^2 \,\widehat{\text{Var}}(\hat{u}_h)}{(\sum_h T_h)^2}\]

strata_df <- data.frame(
  stratum = c("weekday", "weekend"),
  T_h     = c(120L, 80L),
  C_h     = c(280.0, 170.0),
  se_C_h  = c(28.0, 22.0),
  n_h     = c(110L, 70L),
  m_h     = c(9L, 6L)
)

result_strat <- estimate_exploitation_rate(strata = strata_df, by = "stratum")
print(result_strat)
#> 
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: Exploitation Rate (Mark-Recapture)
#> Variance: delta
#> Confidence level: 95%
#> Grouped by: stratum
#> Unit: proportion
#> 
#> # A tibble: 3 × 9
#>   stratum  estimate     se ci_lower ci_upper     n     T     C     m
#>   <chr>       <dbl>  <dbl>    <dbl>    <dbl> <int> <int> <dbl> <int>
#> 1 weekday     0.191 0.0639   0.0643    0.318   110   120   280     9
#> 2 weekend     0.182 0.0749   0.0327    0.332    70    80   170     6
#> 3 .overall    0.187 0.0487   0.0914    0.283   180   200   450    15

Pass aggregate = FALSE to suppress the .overall row.

Reporting rate adjustment

If not all harvested tagged fish are reported, supply the reporting rate \(\lambda \in (0, 1]\):

\[\hat{u}_{adj} = \frac{\hat{u}}{\lambda}\]

result_adj <- estimate_exploitation_rate(
  T              = 200L,
  C              = 450.0,
  se_C           = 42.0,
  n              = 180L,
  m              = 15L,
  reporting_rate = 0.80
)
#> ℹ `C` was supplied as a bare number, so its target cannot be verified.
#>   It must be a season (period) total -- a sampled-day total understates `u` by
#>   the sampling fraction.
#>   Pass the `estimate_total_harvest()` result itself to have this checked.
print(result_adj)
#> 
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: Exploitation Rate (Mark-Recapture)
#> Variance: delta
#> Confidence level: 95%
#> Unit: proportion
#> 
#> # A tibble: 1 × 8
#>   estimate     se ci_lower ci_upper     n     T     C     m
#>      <dbl>  <dbl>    <dbl>    <dbl> <int> <int> <dbl> <int>
#> 1    0.234 0.0619    0.112    0.357   180   200   450    15

Uncertainty in \(\lambda\) is not propagated; treat reporting rate accuracy as a separate sensitivity analysis.