Estimate total harvest by combining effort and HPUE

Computes total harvest estimates by multiplying effort × HPUE with variance propagation via the delta method. Requires a creel design with both count data (for effort estimation) and interview data (for HPUE estimation).

Usage

estimate_total_harvest(
  design,
  by = NULL,
  variance = "taylor",
  conf_level = 0.95,
  target = c("sampled_days", "stratum_total", "period_total"),
  use_trips = NULL,
  estimator = NULL,
  truncate_at = 0.5,
  aggregate_sections = TRUE,
  missing_sections = "warn",
  ci_method = c("delta", "bootstrap"),
  product_variance = c("goodman", "first_order"),
  ci_type = c("symmetric", "log"),
  missing_rate = c("error", "exclude")
)

Arguments

design

A creel_design object with both counts (via add_counts) and interviews (via add_interviews) attached. Both count and interview survey objects must exist. Interview data must include harvest column (specified via harvest parameter in add_interviews).

by

Optional tidy selector for grouping variables. When specified, must match across both effort and HPUE estimates (same calendar strata or interview variables). Accepts bare column names, multiple columns, or tidyselect helpers.

Two kinds of column are not groupings and are refused: the interview id registered by add_catch(), which holds one value per interview and so leaves no within-group variance to estimate, and columns the package derived rather than the user supplying, such as .angler_effort. A wildcard selector drops the derived columns silently; naming one is an error. A column of your own is never treated as derived, whatever it is called.

variance

Character string specifying variance estimation method: "taylor" (default), "bootstrap", or "jackknife". Applied to BOTH effort and HPUE estimation, then combined via delta method.

conf_level

Numeric confidence level (default: 0.95)

target

Character string specifying the effort domain supplied to estimate_effort(). Options are "sampled_days" (default), "stratum_total", or "period_total". This controls which effort domain is multiplied by HPUE so total harvest stays aligned with the requested temporal target. Bus-route and ice designs accept only "sampled_days", as estimate_effort() does: they have no calendar expansion yet (#520).

use_trips

Character. Which interviews contribute to HPUE. "complete" uses only completed trips; "all" includes incomplete ones. Default NULL means "not specified", which resolves to "complete". An interview taken mid-trip reports the harvest so far against the effort so far, and the two do not scale together over the trip, so "all" gives a length-biased rate and a total built from it. Ignored when the design carries no trip status column.

Since GH #271 a roving design routes to all-trip mean-of-ratios here, as it does for estimate_total_catch(), because estimate_harvest_rate() gained the same estimator selection. Both resolve through the same rule, so the total always agrees with its own rate function.

estimator

Character string selecting the rate estimator used for the HPUE component: "ratio-of-means", "mor", or "mortr". Default NULL means "not specified"; see estimate_harvest_rate() for how the pair resolves and when the roving auto-route applies. Bus-route and ice designs accept only "ratio-of-means", because their total is a ratio of Horvitz-Thompson totals with no mean-of-ratios form.

truncate_at

Numeric minimum trip duration in hours for MOR, or NULL to disable truncation. Default 0.5 (30 minutes) per Hoenig et al. (1997). Truncation is not a tuning knob: the untruncated mean-of-ratios estimator has infinite variance. Ignored under ratio-of-means.

aggregate_sections

Logical. When the design was created with add_sections, should a .lake_total row be appended that sums the per-section estimates? Default TRUE. Set to FALSE to return only the per-section rows without the lake total.

missing_sections

Character(1). Action when a registered section is absent from either count data or interview data: "warn" (default) inserts an NA row with data_available = FALSE, "error" raises a hard error.

ci_method

character. "delta" (default) returns only delta-method CIs. "bootstrap" additionally returns ci_lo_boot/ci_hi_boot using survey bootstrap resampling. Only applies to bus-route/ice designs.

product_variance

character. Variance formula for the product $E \times H$. "goodman" (default) uses Goodman's (1960) unbiased estimator $E^2 Var(H) + H^2 Var(E) - Var(E)Var(H)$; "first_order" omits the cross-term, which is conservative. Both assume $E$ and $H$ are independently estimated. When both components are so imprecise that the subtraction would give a non-positive variance, the first-order value is used as a floor.

ci_type

character. Shape of the confidence interval. "symmetric" (default) gives $\hat\theta \pm z \cdot SE$ clamped at zero. "log" applies a log-transform for a strictly positive CI.

missing_rate

character. What to do when a stratum (or stratum x by cell) has effort but no rate – no interviews, or only trips the rate cannot use. Its harvest is unknown, not zero. "error" (default) stops and names the cells and their share of effort. "exclude" reports the total over the covered cells only, warns, gives a by group with no covered cell an NA row, and records what it left out in excluded_strata on the result. Not available for bus-route and ice designs, whose totals have no per-stratum product (GH #373). A cell whose effort is a known zero (estimate and SE both 0) is not counted as missing: its product is 0 for any rate. An ungrouped total on a sectioned design pools one rate across each section's strata, so an uncovered stratum there is refused under either setting; use by on the stratum column to get the covered cells and the "exclude" record.

Value

A creel_estimates S3 object with method = "product-total-harvest". The estimator component records the rate estimator this total is a product of, as you asked for it: method names the product form and is the same string whichever estimator produced it. For bus-route and ice designs, returns a bus-route HT estimate with method = "ht-total-harvest" and a "site_contributions" attribute.

For sectioned designs the per-section rows carry prop_of_lake_total, the section's share of the lake-wide total, and se_prop_of_lake_total, its standard error. The share is a ratio whose numerator is one of its own denominator's terms, and whose numerator and denominator are each products of an effort and a rate estimated from different designs, so the error is derived by delta method from the same section variances and covariance the .lake_total row's own standard error is built from. The .lake_total row reports se_prop_of_lake_total = 0: its share of itself is exactly 1 by construction and was never estimated. A section with no data reports NA for both. Neither column is produced on the grouped path.

Details

Total harvest is computed as Effort × HPUE. Variance is propagated using the delta method, which accounts for uncertainty in both estimates. The formula for independent estimates is approximately:

$$Var(E \times H) \approx E^2 \cdot Var(H) + H^2 \cdot Var(E)$$

Variance is computed via a stratified delta-method sum in compute_stratum_product_sum(), not via survey::svycontrast().

Sectioned designs: When add_sections has been called on the design, each section is estimated independently. The lake-wide total is sum(TH_i), not E_total * HPUE_pooled. The lake-wide SE adds, to the sections' own variances, the covariance of their effort: sections counted on the same days share a day effect, so their effort totals move together. That term is 2 * sum_h R_ih R_jh Cov(E_ih, E_jh) over pairs of sections and strata, estimated from the days both sections were counted; sections never counted on the same day add nothing (GH #414). Rates come from each section's own interviews and are treated as independent across sections. Sections declared with add_sections(shared_count_times = TRUE) also add the within-day covariance of their counts, the term the effort lake row carries for them.

by = <species> is supported on a sectioned design: harvest is apportioned against each section's own whole effort, giving one row per section per species. As with any other grouping, the sectioned result then carries no .lake_total row and no prop_of_lake_total.

Design compatibility requirements:

Count data must be attached via add_counts() for effort estimation
Interview data must be attached via add_interviews() for HPUE estimation
Harvest column must be specified in add_interviews (harvest parameter)
Grouped estimation requires identical grouping variables for both estimates
Calendar stratification must be shared between counts and interviews

Why there is no `targeted` argument

The rate functions accept targeted = FALSE, which restricts the domain to the interviews that recorded some of the species being estimated. The totals deliberately do not, because a total is a rate multiplied by an effort base and the two would no longer describe the same set of trips.

A targeted rate is conditional on having recorded the species; total effort is not. Multiplying one by the other applies a conditional rate to an unconditional base. On the package's own example data one species' rate is 0.48 fish/hr over all 50 trips and 2.00 fish/hr over the 12 that caught it, so expanding the targeted rate by total effort returns roughly 223 fish where 30 were actually caught.

The domain-consistent product — the targeted rate times the effort of the trips that recorded the species — is well defined in the sample but cannot be expanded: it needs the season-wide effort of species-catching trips, which no creel design observes.

So a targeted rate is available and a targeted total is not, and that is a property of the estimand rather than a gap in the implementation (GH #307).

What the pooled total assumes

Effort comes from the counts, so a total can only be broken down by an attribute the counts classify. When a domain appears in the interviews but not in the counts, the only available total is E_total * rate_pooled, where the pooled rate is a ratio of means weighted by the interview sample's composition over that domain. Had the domain been classified in the counts it would be a stratum and the total would be sum(E_h * rate_h), which is unbiased whatever the interview composition happens to be.

The two agree only when the interview sample's effort composition matches the true effort composition, and interview selection is not proportional to effort by construction of the standard designs. Access interviews intercept completed trips, over-representing anglers who must return to a fixed point: Malvestuto (1996) notes that it is “usually impossible to sample all angler types proportional to their level of effort”, a particular problem for bank anglers who may be “widely dispersed along the shoreline and not associated with well-defined access sites”. Roving interviews are length-biased toward longer trips. So the mix differs by design rather than by accident, and where levels differ in rate the pooled total inherits that difference.

None of this is verifiable from within the data, because the counts carry no composition to compare against. Where it is detectable – the interviews hold an unclassified categorical domain and the crude rate differs materially across its levels – a warning of class creel_warning_pooled_domain_mix is raised. It flags a risk, not a defect. Classifying the domain in the count data is what removes the assumption.

Unit of the total

The reported unit is derived from the two factors, never declared. A total is "fish" only when a per-angler-hour rate multiplies an effort in angler-hours; anything else reports NA_character_, meaning unknown.

Two ways to fail to cancel:

The effort unit is unknown. design$effort_unit is NA whenever add_counts() received no period_length_col, because a bare count column may be an instantaneous head count or effort the caller already expanded, and nothing can tell the two apart. Unknown times known is unknown. Supply period_length_col to make the total's unit derivable. An aerial design is the exception: it refuses period_length_col, but h_open is its period length, so its effort is recorded as angler-hours and its totals are "fish" (GH #478, #474).
The denominators disagree. A rate per party-hour times an effort in angler-hours is not a count of fish. Pass n_anglers to add_interviews() so the rate is per angler-hour.

The estimate itself is unaffected in both cases – only the label changes. Until version 5.2.0 the unit was the literal "fish" regardless of either factor (GH #213).

Examples

library(tidycreel)
data(example_calendar)
data(example_counts)
data(example_interviews)

# Create design with both counts and interviews including harvest
design <- creel_design(example_calendar, date = date, strata = day_type)
design <- add_counts(design, example_counts)
design <- add_interviews(design, example_interviews,
  catch = catch_total, harvest = catch_kept, effort = hours_fished,
  n_anglers = n_anglers,
  trip_status = trip_status, trip_duration = trip_duration
)
#> ℹ Added 22 interviews: 17 complete (77%), 5 incomplete (23%)

# Estimate total harvest
total_harvest <- estimate_total_harvest(design)
#> Warning: Instantaneous counts were expanded without a period length.
#> ℹ No `period_length_col` was supplied to `add_counts()`, so the estimate is the
#>   count column summed over days.
#> ! If that column holds an instantaneous angler count, the result is in
#>   angler-days, not angler-hours.
#> ℹ Supply the period each count was randomised within: `add_counts(design,
#>   counts, period_length_col = <col>)`.
#> ℹ Repeats on this design are recorded, not printed: see `get_warnings()`.
#> Warning: `estimate_total_harvest()` is pooling over a domain the counts do not classify:
#> angler_method.
#> ! The rate differs across its levels in these interviews (angler_method:
#>   artificial 0.667, bait 0.483, fly 0.897), so the total depends on the
#>   interview sample's mix over that domain.
#> ℹ Without the domain in the counts the total is `E_total * rate_pooled`,
#>   weighted by the interview mix rather than the effort mix. Interview selection
#>   is not proportional to effort by construction (Malvestuto 1996).
#> ℹ This is a risk, not an error: the counts carry no composition to check
#>   against, so it cannot be verified from the data.
#> ℹ Classifying angler_method in the count data removes the assumption -- the
#>   total becomes `sum(E_h * rate_h)`.
#> ℹ Repeats on this design are recorded, not printed: see `get_warnings()`.
#> Warning: 1 rate cell in "estimate_total_harvest" rests on fewer than 10 trips.
#> • day_type=weekend: n=7
#> ℹ The total multiplies each cell's effort by a rate built from that cell's
#>   trips alone; a ratio from so few trips is unstable, and the rate functions
#>   refuse a group this small.
#> ℹ Combine cells or collect more interviews. Totals warn here rather than refuse
#>   (GH #417).
#> ℹ Repeats on this design are recorded, not printed: see `get_warnings()`.
print(total_harvest)
#> 
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: Total Harvest (Effort × HPUE)
#> Variance: Taylor linearization
#> Confidence level: 95%
#> Effort target: sampled_days
#> 
#> # A tibble: 1 × 5
#>   estimate    se ci_lower ci_upper     n
#>      <dbl> <dbl>    <dbl>    <dbl> <int>
#> 1     225.  28.1     165.     284.    17

# Compare components
effort_est <- estimate_effort(design)
hpue_est <- estimate_harvest_rate(design)
#> ℹ Filtering to complete trips for HPUE estimation
#>   (n=17, 77.3% of 22 interviews) [default]
#> Warning: Small sample size for harvest estimation.
#> ! Sample size is 17. Ratio estimates are more stable with n >= 30.
#> ℹ Variance estimates may be unstable with n < 30.
#> ℹ Repeats on this design are recorded, not printed: see `get_warnings()`.
# The total is close to, but not exactly, effort times HPUE
c(
  total = total_harvest$estimates$estimate,
  effort_x_hpue = effort_est$estimates$estimate * hpue_est$estimates$estimate
)
#>         total effort_x_hpue 
#>      224.5239      221.7674 

# Grouped estimation needs at least 10 interviews per group, so check
# the sample sizes before grouping
table(design$interviews$day_type)
#> 
#> weekday weekend 
#>      13       9