This page maps each formula in tidycreel’s bus-route estimators to its published source and to the function that computes it. It is a reference for checking the implementation, not a tutorial; the bus-route How-to page shows the workflow.
The sources, and what each supplies:
Code references below name functions, not line numbers, which go stale.
Formula:
\[\pi_i = p\_site_i \times p\_period\]
p_period is the probability that the sampled period (for
example the morning or the afternoon) is the one worked; it is constant
within a circuit. p_site is the site’s weight within the
circuit, given in one of two ways:
| Argument(s) | p_site |
Constraint | Design |
|---|---|---|---|
wait_time, circuit_time |
\(w_i / T\) | waits add up to no more than \(T\) | bus route: every site visited on each pass |
p_site |
as given | adds up to 1 within a circuit | a randomly chosen site on each pass |
For a bus route, \(w_i\) is the scheduled wait at site \(i\) and \(T\) “the total time to complete a full circuit of the route, including traveling and waiting” (Pollock et al. 1994, p. 238). A party at the site can be seen only while the crew is there, so \(w_i/T\) is the share of the period during which it can be observed. The shares add up to less than 1 whenever the route has travel time. Rescaling them to add up to 1 divides by the total wait instead of \(T\), and lowers every total by the travelling share of the circuit (GH #539).
R implementation: creel_design()
computes .p_site (from waits) and .pi_i in the
sampling frame and checks the constraints; add_interviews()
joins .pi_i to each interview by site and circuit.
Formula:
\[\text{expansion}_i = \frac{n\_counted_i}{n\_interviewed_i}\]
When not every party at a site is interviewed, the interviewed parties stand for those that were not. With every party interviewed, the expansion is 1, as in the examples below.
R implementation: add_interviews()
stores .expansion. Two edge cases:
estimate_effort() sets
it to 0 deliberately.add_interviews() warns. Their effort and
catch are unknown, not zero, so the effort, catch, harvest and release
totals, and the catch, harvest and release rates, stop with an error
that names each such site-day. Remove those rows to leave the parties
out deliberately, or impute their values and set
n_interviewed to the number of parties the values stand
for. The same refusal covers a missing count or a missing value on an
interviewed row.Formula (Jones & Pollock 2012, Eq. 19.4):
\[\hat{E} = \sum_{i} \frac{e_i}{\pi_i}, \qquad e_i = a_i \times \text{expansion}_i\]
where \(a_i\) is interview \(i\)’s angler effort: the hours fished times
n_anglers, or the hours fished alone when no party size is
given. Pollock et al. write \(\pi_j\)
for the probability of the sampled period, tidycreel’s
p_period. With \(\pi_i =
(w_i/T)\,\pi_j\) the sum is their Eq. 15.10 term for term:
\[\hat{E} = T \sum_i \frac{1}{w_i} \sum_j \frac{e_{ji}}{\pi_j}\]
Unit: angler-hours when
add_interviews() was given n_anglers;
party-hours when it was not (it warns). The effort column is read as
hours.
R implementation: estimate_effort()
dispatches to estimate_effort_br(), which forms
.e_i and .contribution \(= e_i/\pi_i\) and sums them.
get_site_contributions() returns the per-interview
table.
Formula (Jones & Pollock 2012, Eq. 19.5; Pollock et al. 1994, Eq. 15.11):
\[\hat{H} = \sum_i \frac{h_i}{\pi_i}, \qquad h_i = \text{harvest}_i \times \text{expansion}_i\]
and likewise for catch and released fish. These totals are Horvitz–Thompson totals of the fish themselves. They are not effort times a rate, so they need no assumption about the covariance of the two.
R implementation:
estimate_total_catch(),
estimate_total_harvest() and
estimate_total_release() dispatch to
estimate_total_catch_br(),
estimate_total_harvest_br() and
estimate_total_release_br(), which share
br_build_estimates().
Eq. 19.4 and 19.5 are totals; the sources define no separate bus-route rate. The rate this design supports is the ratio of the two totals:
\[\widehat{\text{HPUE}} = \frac{\hat{H}}{\hat{E}} = \frac{\sum_i h_i/\pi_i}{\sum_i e_i/\pi_i}\]
This is a ratio of means, reported as
method = "ratio-of-means-hpue" (or -cpue,
-rpue). Its unit follows the denominator: fish per
angler-hour, or per party-hour without n_anglers.
Variance. survey::svyratio() linearises
numerator and denominator together. The two totals come from the same
interviews and are strongly positively correlated, so dividing two
separately estimated totals and treating them as independent overstates
the SE.
R implementation:
estimate_harvest_rate(), estimate_catch_rate()
and estimate_release_rate() dispatch to
estimate_harvest_br() (catch and release reuse it with
their own column), whose ratio step is
br_harvest_rate_estimates(). On a bus-route design
use_trips defaults to "complete".
Anglers interviewed mid-trip need a different estimator, selected
with use_trips = "incomplete".
Hoenig, Jones, Pollock, Robson & Wade (1997) compare the two candidates under roving-type sampling: encounter probability proportional to time fished, and catch observed only up to the interview. Writing \(C_j\) for the catch and \(L_j\) for the hours fished by angler \(j\) at the time of interview, the ratio of means has expectation (their Eq. 3.2)
\[E[R_1] \;\approx\; \frac{\sum_j \lambda_j L_j^{*2}}{\sum_j L_j^{*2}}\]
a weighting of individual rates \(\lambda_j\) by the square of completed trip length \(L_j^{*}\). So \(R_1\) “does not provide an estimate of catch rate that can be used with an independent estimate of total effort to provide an unbiased estimate of total catch except in the unrealistic case where \(\lambda\) is constant over all anglers.” The mean of ratios has approximately the correct expectation.
On a bus route the interviews are not equally likely, so each enters a Hájek weighted mean:
\[\widehat{\text{HPUE}}_{\text{inc}} \;=\; \frac{\sum_i w_i \, (h_i / a_i)}{\sum_i w_i}, \qquad w_i = \frac{\text{expansion}_i}{\pi_i} \times n\_anglers_i\]
where \(a_i\) is angler-effort. The
\(n\_anglers_i\) factor applies when
party size varies: an interview is a party, the rate is per angler, and
meeting a party meets all its anglers (GH #513). It is reported as
method = "mean-of-ratios-hpue".
Truncation. The mean of ratios has infinite
asymptotic variance as trip length goes to zero. Hoenig et al. recommend
disregarding trips shorter than 30 minutes, which is the
truncate_at = 0.5 default. The threshold applies to elapsed
trip duration, not to angler-hours: a party of five fishing twelve
minutes supplies a full angler-hour and is still the unstable case.
truncate_at = NULL disables truncation and warns.
use_trips = "diagnostic" returns both rates. They are
different estimators, each the one its trip type supports, but both
divide by the same effort column, so they share a unit (per angler-hour
with n_anglers, per party-hour without) and the gap between
them reflects trip status rather than a change of unit.
R implementation:
br_incomplete_harvest_rate(), called from
estimate_harvest_br().
Design. The per-interview contributions are analysed as a stratified sample with the sampled day as the primary sampling unit, stratified by the design’s strata. Interviews on the same day share that day’s conditions, so they are not independent draws (GH #198).
Method. survey::svytotal() for totals
and survey::svyratio() for rates, by Taylor linearisation
(Lumley 2010); variance = "bootstrap" or
"jackknife" replicate the same design. A stratum with one
sampled day has no between-day variance and is refused.
Within-day variance. Robson & Jones (1989, p. 95) show that with one crew per route per day the within-day component cannot be estimated separately. The day-level variance contains it, so the SE includes it, but tidycreel reports no separate within-day term for bus routes.
Target. The total covers the sampled
days. Robson & Jones (Eq. 13) expand the mean daily
estimate to all \(N\) days of a
stratum; tidycreel does not implement that expansion for bus routes, and
target = "stratum_total" or "period_total" is
refused rather than ignored (GH #520). The result records
effort_target = "sampled_days".
R implementation: br_add_psu_key() sets
the PSU; get_variance_design() applies the variance
method.
Pollock et al. give one day of direct-expansion interviews on a five-site route: a 360-minute circuit with waits of 30, 30, 15, 60 and 60 minutes, and the morning or afternoon chosen with probability 0.5. Site 3 had no anglers. They report \(\hat{E} = 24{,}814.80\) minutes (Table 15.16) and \(\hat{C} = 227.88\) fish (Table 15.17).
A design needs two sampled days for a between-day variance, so the published day is entered twice. Each day’s estimate is then the published one, and the sampled-day total is twice it. Trip durations are converted to hours.
library(tidycreel)
frame <- data.frame(
site = as.character(1:5),
circuit = "am_pm",
wait_min = c(30, 30, 15, 60, 60),
circuit_min = 360
)
one_day <- data.frame(
site = c("1", "1", "2", "2", "4", "4", "4", "5"),
trip_min = c(135, 165, 175, 215, 150, 145, 168, 225),
catch = c(3, 2, 0, 1, 1, 2, 3, 1)
)
days <- as.Date(c("2024-06-01", "2024-06-02"))
interviews <- do.call(rbind, lapply(days, function(d) cbind(date = d, one_day)))
interviews$day_type <- "weekday"
interviews$circuit <- "am_pm"
interviews$trip_hours <- interviews$trip_min / 60
interviews$trip_status <- "complete"
# Every party at a site was interviewed.
interviews$n_counted <- ave(rep(1L, nrow(interviews)), interviews$date, interviews$site,
FUN = length
)
interviews$n_interviewed <- interviews$n_countedPollock et al.’s estimate is party effort, the trip length of each
party, so the interviews are added without n_anglers.
add_interviews() warns that it is assuming one angler per
interview, and the estimates are labelled party-hours,
which is what this example is.
design <- creel_design(
data.frame(date = days, day_type = "weekday"),
date = date, strata = day_type, survey_type = "bus_route",
sampling_frame = frame, site = site, circuit = circuit,
wait_time = wait_min, circuit_time = circuit_min, p_period = 0.5
)
design <- add_interviews(design, interviews,
effort = trip_hours, catch = catch, trip_status = trip_status,
n_counted = n_counted, n_interviewed = n_interviewed
)
#> Warning: ! No `n_anglers` provided — assuming 1 angler per interview.
#> ℹ Pass `n_anglers = <column>` to use actual party sizes for angler-hour
#> normalization.
#> ℹ If the interviews really are one angler each, pass `n_anglers = 1` to state
#> that and silence this warning.
#> ℹ 2 interviews have zero catch.
#> Zero catch is usual; check it is not missing data recorded as 0.
#> ℹ Added 16 interviews: 16 complete (100%), 0 incomplete (0%)effort <- estimate_effort(design)
effort$estimates
#> # A tibble: 1 × 5
#> estimate se ci_lower ci_upper n
#> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 827. 0 827. 827. 16
effort$unit
#> [1] "party-hours"
# Per day, in minutes: Pollock et al. print 24,814.80 from site terms rounded to two
# decimals; unrounded it is 24,816.
effort$estimates$estimate / 2 * 60
#> [1] 24816
catch <- estimate_total_catch(design)
# Per day: Pollock et al. print 227.88 from the rounded 0.633; unrounded it is 228.
catch$estimates$estimate / 2
#> [1] 228The standard error is 0 only because the two entered days are
identical, so there is no between-day variation to measure. It says
nothing about Pollock et al.’s survey, which reports one day and no
variance. tests/testthat/test-bus-route-wait-time.R checks
both figures.
| Quantity | Formula | Source | R function |
|---|---|---|---|
| Site weight (bus route) | \(p\_site_i = w_i / T\) | Robson & Jones (1989) Eq. 12; Pollock et al. (1994) Eq. 15.10 (both as \(T/w_i\)) | creel_design(wait_time, circuit_time) |
| Inclusion probability | \(\pi_i = p\_site_i \times p\_period\) | Pollock et al. (1994) Eq. 15.10 (\(\pi_j\)) | creel_design(), add_interviews() |
| Enumeration expansion | \(n\_counted / n\_interviewed\) | — | add_interviews() |
| Effort total | \(\hat{E} = \sum e_i/\pi_i\) | Jones & Pollock (2012) Eq. 19.4 | estimate_effort_br() |
| Catch/harvest/release total | \(\hat{H} = \sum h_i/\pi_i\) | Jones & Pollock (2012) Eq. 19.5; Pollock et al. (1994) Eq. 15.11 | br_build_estimates() |
| Rate, complete trips | \(\hat{H}/\hat{E}\) | ratio of Eq. 19.5 to Eq. 19.4 | br_harvest_rate_estimates() |
| Rate, incomplete trips | Hájek mean of \(h_i/a_i\) | Hoenig et al. (1997) | br_incomplete_harvest_rate() |
| Variance | day-PSU stratified design, Taylor or replicates | Lumley (2010) | br_add_psu_key(),
get_variance_design() |
Each site’s contribution is divided by its \(\pi_i\), so an error in \(\pi_i\) is an error in that site’s share of the total. Two common errors, on the Pollock et al. route above:
Hoenig, J. M., Jones, C. M., Pollock, K. H., Robson, D. S., & Wade, D. L. (1997). Calculation of catch rate and total catch in roving surveys of anglers. Biometrics, 53(1), 306–317.
Jones, C. M., & Pollock, K. H. (2012). Recreational angler survey methods: estimation of effort, harvest, and released catch. Chapter 19 in Fisheries Techniques (3rd ed.), pp. 883–919. American Fisheries Society.
Lumley, T. (2010). Complex Surveys: A Guide to Analysis Using R. Wiley.
Pollock, K. H., Jones, C. M., & Brown, T. L. (1994). Angler survey methods and their applications in fisheries management. American Fisheries Society Special Publication 25.
Robson, D., & Jones, C. M. (1989). The theoretical basis of an access site angler survey design. Biometrics, 45, 83–98.