---
title: "Bus-Route Survey: Equation Traceability"
output:
  rmarkdown::html_vignette:
    highlight: null
vignette: >
  %\VignetteIndexEntry{Bus-Route Survey: Equation Traceability}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Overview

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](bus-route-surveys.html) shows the workflow.

The sources, and what each supplies:

- **Robson & Jones (1989)** define the bus-route design and its effort estimator: the
  crew visits every site on each pass, waits a scheduled time at each, and expands what it
  sees by the circuit time over the wait (their Eq. 12).
- **Pollock, Jones & Brown (1994, section 15.5.3)** give the same estimator as Eq. 15.10
  (effort) and 15.11 (catch), with worked examples in Tables 15.14–15.17.
- **Jones & Pollock (2012, Eq. 19.4–19.5, pp. 911–912)** give the general
  Horvitz–Thompson totals, $\hat{E} = \sum e_i/\pi_i$ and $\hat{H} = \sum h_i/\pi_i$,
  "where $\pi_i$ is the total probability that the ith sampling unit is included in the
  sample". tidycreel writes the bus-route estimator in this form.
- **Hoenig et al. (1997)** supply the estimator for anglers interviewed mid-trip.

Code references below name functions, not line numbers, which go stale.

---

## 1. Site Weight and Inclusion Probability

**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.

---

## 2. Enumeration Expansion

**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:

- **Nobody seen** ($0/0$): the site-visit contributes zero effort, because no party was
  there; `estimate_effort()` sets it to 0 deliberately.
- **Parties seen, none interviewed** ($n\_counted > 0$, $n\_interviewed = 0$): the
  expansion is undefined and `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.

---

## 3. Effort

**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.

---

## 4. Catch, Harvest and Release Totals

**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()`.

---

## 5. Rates from Complete Trips

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"`.

---

## 6. Rates from Incomplete Trips

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()`.

---

## 7. Variance and What the Total Covers

**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.

---

## 8. Worked Check: Pollock et al. (1994), Tables 15.15–15.17

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.

```{r pollock-frame}
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_counted
```

Pollock 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.

```{r pollock-design}
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
)
```

```{r pollock-estimates}
effort <- estimate_effort(design)
effort$estimates
effort$unit

# 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

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
```

The 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.

---

## 9. Summary Traceability Table

| 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()` |

---

## 10. Why the Weights Matter

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:

- **A fixed $\pi_i = 0.5$ for every site.** That is not an inclusion probability of the
  design. Each site is then off by the ratio of 0.5 to its true $\pi_i$; site 1's true
  $\pi_i$ is $30/360 \times 0.5 = 0.042$, so it would be understated twelvefold.
- **Wait shares rescaled to add up to 1** ($w_i / \sum w$). Every total is multiplied by
  $\sum w / T = 195/360 = 0.54$: the 165 minutes of travel are treated as if they were
  spent watching.

---

## References

- 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.