This page explains how tidycreel
turns interviews into a catch rate, why the rate estimator depends on
how the interviews were collected, and how a total catch and its
standard error are built from effort and rate. Effort itself is
explained in vignette("effort-pipeline").
interviews -> catch rate per stratum -> total = sum over strata of (effort x rate)
(access: ratio of means; (standard error: Goodman's product
roving: mean of ratios) variance, added over strata)Total catch is effort times catch per unit of effort. For that product to be the total, the rate has to be total catch divided by total effort over all anglers in the stratum, in fish per angler-hour. That is the quantity every estimator below is judged against. It is not the same as how well the average angler did: an angler who fished one hour counts as much as one who fished eight in an average of anglers, but contributes an eighth as much effort.
From \(n\) interviews with catch \(c_i\) and hours fished \(h_i\):
\[\hat{R}_{\text{ROM}} = \frac{\sum c_i}{\sum h_i} \qquad\qquad \hat{R}_{\text{MOR}} = \frac{1}{n} \sum \frac{c_i}{h_i}\]
The ratio of means (ROM) is total catch over total hours. The mean of ratios (MOR) is the average of each party’s own rate. On the same ten interviews they differ:
catch <- c(2, 0, 6, 1, 3, 4, 0, 2, 1, 5)
hours <- c(4, 1, 3, 2, 2, 4, 2, 2, 1, 3)
sum(catch) / sum(hours) # ratio of means
#> [1] 1
mean(catch / hours) # mean of ratios
#> [1] 0.9166667Which one is right depends on how the interviewed anglers were selected, because that decides what each one estimates (Pollock et al. 1997).
At an access point the clerk interviews parties as they finish, so every trip has the same chance of being interviewed whatever its length. The ratio of means then estimates total catch over total effort, the quantity a total needs. The mean of ratios estimates the average angler’s rate instead, a different quantity (Pollock et al. 1997).
A roving clerk meets parties while they fish. A party on the water for eight hours is eight times as likely to be met as one fishing for one hour, so the interviews over-represent long trips (length-biased sampling). Under that selection the mean of ratios, computed on the catch and hours so far, estimates total catch over total effort approximately, and the ratio of means does not: it gives the long trips extra weight a second time (Pollock et al. 1997; Hoenig et al. 1997).
That argument is about anglers, and an interview is a party. Meeting
a party meets every angler in it, so a party of three stands for three
anglers, each with the party’s rate per angler-hour. With party sizes
\(m_i\) (n_anglers) and
the rate per angler-hour \(r_i = c_i / (h_i
m_i)\), the mean of ratios is
\[\hat{R}_{\text{MOR}} = \frac{\sum m_i r_i}{\sum m_i}\]
Averaging the parties with equal weight would estimate a different
quantity whenever the rate per angler differs between small and large
parties. Without n_anglers the rate is per party-hour and
each party counts once, which is the formula above.
Parties met a few minutes into a trip have rates that are mostly
noise (one fish in ten minutes is six an hour), which inflates the
variance of the mean of ratios. Leaving out trips shorter than a minimum
lowers the variance with negligible bias; Pollock et al. (1997) suggest
20–30 minutes, and tidycreel’s default is 30
(truncate_at = 0.5).
The estimator follows the interview type declared on the design:
interview_type |
Default estimator | Trips used |
|---|---|---|
"access" (default) |
ratio of means | complete trips |
"roving" |
mean of ratios, trips under 0.5 h left out | all trips |
estimator and use_trips override the
default. vignette("instantaneous-counts") runs both
interview types on one season; vignette("incomplete-trips")
covers using incomplete trips from an access survey.
The season from vignette("tidycreel"): 120 days in two
strata, instantaneous counts, and access interviews.
design <- creel_design(example_season_calendar, date = date, strata = day_type)
design <- add_counts(design, example_season_counts,
count_col = total_anglers,
period_length_col = daylight_hours
)
design <- add_interviews(design, example_season_interviews,
catch = catch_total,
effort = hours_fished,
trip_status = trip_status,
n_anglers = n_anglers
)
effort <- estimate_effort(design, by = day_type, target = "period_total")
rate <- estimate_catch_rate(design, by = day_type)
effort
#>
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: Total
#> Variance: Taylor linearization
#> Confidence level: 95%
#> Grouped by: day_type
#> Effort target: period_total
#> Unit: angler-hours
#>
#> # A tibble: 2 × 8
#> day_type estimate se se_between se_within ci_lower ci_upper n
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 weekday 7912 1063. 1063. 0 5753. 10071. 28
#> 2 weekend 2082. 725. 725. 0 609. 3556. 8
rate
#>
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: Ratio-of-Means CPUE
#> Variance: Taylor linearization
#> Confidence level: 95%
#> Grouped by: day_type
#> Unit: fish/angler-hour
#>
#> # A tibble: 2 × 6
#> day_type estimate se ci_lower ci_upper n
#> <chr> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 weekday 0.208 0.0453 0.119 0.297 138
#> 2 weekend 0.280 0.102 0.0794 0.481 34Weekend parties caught fish faster than weekday parties, and weekdays carry most of the effort. The total is built within each stratum and then added:
\[\hat{C} = \sum_h \hat{E}_h \hat{R}_h\]
total <- estimate_total_catch(design, target = "period_total")
total
#>
#> ── Creel Survey Estimates ──────────────────────────────────────────────────────
#> Method: Total Catch (Effort × CPUE)
#> Variance: Taylor linearization
#> Confidence level: 95%
#> Effort target: period_total
#> Unit: fish
#>
#> # A tibble: 1 × 5
#> estimate se ci_lower ci_upper n
#> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 2230. 506. 1231. 3230. 172
E <- effort$estimates
R <- rate$estimates
sum(E$estimate * R$estimate)
#> [1] 2230.484The by-hand sum matches. Multiplying the season effort by a single season rate would give a different number, because the pooled rate weights each stratum by its interviews, not by its effort:
estimate_effort(design, target = "period_total")$estimates$estimate *
estimate_catch_rate(design)$estimates$estimate
#> [1] 2205.907Here the difference is about 1%, but it grows when strata differ more in rate and in how heavily they were interviewed.
Within a stratum, the total is a product of two estimates, each with its own error. If the two are independent, the exact variance of a product is (Goodman 1960, equation 2)
\[\text{Var}(\hat{E}\hat{R}) = E^2\,\text{Var}(\hat{R}) + R^2\,\text{Var}(\hat{E}) + \text{Var}(\hat{E})\,\text{Var}(\hat{R})\]
The unknown \(E\) and \(R\) have to be replaced by their estimates. Because \(\hat{E}^2\) overstates \(E^2\) by \(\text{Var}(\hat{E})\) on average, and likewise for \(\hat{R}^2\), each of the first two terms is biased upward by \(\text{Var}(\hat{E})\,\text{Var}(\hat{R})\), together overshooting the exact variance by that amount. Goodman’s unbiased estimator (equation 5) therefore subtracts the third term instead of adding it:
\[\widehat{\text{Var}}(\hat{E}_h\hat{R}_h) = \hat{E}_h^2\,\widehat{\text{Var}}(\hat{R}_h) + \hat{R}_h^2\,\widehat{\text{Var}}(\hat{E}_h) - \widehat{\text{Var}}(\hat{E}_h)\,\widehat{\text{Var}}(\hat{R}_h)\]
This is tidycreel’s default
(product_variance = "goodman"). The strata are sampled
independently, so their variances add:
v_goodman <- with(
data.frame(e = E$estimate, se_e = E$se, r = R$estimate, se_r = R$se),
sum(e^2 * se_r^2 + r^2 * se_e^2 - se_e^2 * se_r^2)
)
sqrt(v_goodman)
#> [1] 506.2577which is the se printed for the total. The familiar
two-term delta method leaves the subtraction out; it is available as
product_variance = "first_order" and is slightly larger, so
conservative:
estimate_total_catch(design,
target = "period_total", product_variance = "first_order"
)$estimates$se
#> [1] 513.9317When both estimates are very imprecise the subtraction can make the variance negative; tidycreel then falls back to the first-order value.
Both formulas assume effort and rate are estimated independently, so their covariance is zero. That is an assumption, not a property of the design: counts and interviews are collected on the same sampled days. If busy days were also days of fast fishing, the two estimates would tend to err in the same direction and the standard error would understate the uncertainty; if the reverse, it would overstate it. tidycreel does not estimate this covariance.
Each stratum’s variance has an effort part, \(\hat{R}_h^2 \widehat{\text{Var}}(\hat{E}_h)\), and a rate part, \(\hat{E}_h^2 \widehat{\text{Var}}(\hat{R}_h)\):
data.frame(
day_type = E$day_type,
from_effort = R$estimate^2 * E$se^2,
from_rate = E$estimate^2 * R$se^2
)
#> day_type from_effort from_rate
#> 1 weekday 48936.77 128469.82
#> 2 weekend 41230.16 45489.06In both strata the rate part is the larger: about two and a half times the effort part on weekdays, and only slightly larger on weekends. This says where the uncertainty comes from. It does not by itself say whether to collect more counts or more interviews: that also depends on what each costs, and on the fact that more sampled days improve both.