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

## Overview

This document maps every formula used in tidycreel's bus-route estimation to its
published primary source, including the specific equation number and page number. It is
intended as an auditable record for reviewers, developers, and users who want to verify
that the implementation matches the statistical theory.

The bus-route estimation framework in tidycreel follows two primary sources:

- **Jones & Pollock (2012)**: Chapter 19 of *Fisheries Techniques* (3rd ed.), pp. 883–919.
  This chapter provides the general nonuniform probability creel survey estimators.
- **Malvestuto (1996)**: Chapter 20 of *Fisheries Techniques* (2nd ed.), pp. 591–623.
  Box 20.6 provides a worked numerical example that tidycreel reproduces exactly.

The four key computations are: inclusion probability (πᵢ), enumeration expansion,
expanded effort per interview, and the Horvitz-Thompson total estimators for effort and
harvest. Each is documented in a section below.

---

## 1. Inclusion Probability (πᵢ)

**Published source:** Jones & Pollock (2012), p. 912

**Formula:**

$$\pi_i = p\_site_i \times p\_period$$

where $p\_site_i$ is the probability of selecting site $i$ during a given circuit pass,
and $p\_period$ is the probability that the sampling period is included in the survey.
For uniform period sampling (all periods equally likely), $p\_period$ is the same for all
sites and circuits.

**Statistical meaning:** $\pi_i$ is the *inclusion probability* — the marginal probability
that unit $i$ is included in the sample under the two-stage design (stage 1: select a
period; stage 2: traverse the circuit). Jones & Pollock (2012) treat this as the product
of the two independent selection probabilities.

**R implementation:**

- `creel_design()` in `R/creel-design.R`: accepts `p_site` and `p_period` columns from the
  `sampling_frame` argument and precomputes `pi_i = p_site * p_period` for each
  site×circuit combination. The result is stored in `design$bus_route$sampling_frame`.
- `add_interviews()` in `R/creel-design.R`: joins the precomputed `pi_i` value from the
  sampling frame to each interview row, storing it as the `.pi_i` column.

**Malvestuto (1996) Box 20.6 values:**

| Site | p_site | p_period | πᵢ    |
|------|--------|----------|-------|
| A    | 0.30   | 0.50     | 0.150 |
| B    | 0.25   | 0.50     | 0.125 |
| C    | 0.40   | 0.50     | 0.200 |
| D    | 0.05   | 0.50     | 0.025 |

---

## 2. Enumeration Expansion

**Published source:** Malvestuto (1996), Box 20.6, p. 614

**Formula:**

$$\text{expansion}_i = \frac{n\_counted_i}{n\_interviewed_i}$$

where $n\_counted_i$ is the total number of angler parties observed at site $i$ during
the visit, and $n\_interviewed_i$ is the number of those parties actually interviewed.

**Statistical meaning:** When not all parties at a site are interviewed, the expansion
factor rescales the interviewed sample to represent the full party count. If all parties
are interviewed ($n\_counted = n\_interviewed$), the expansion is 1 and no adjustment is
needed. This is the case in Box 20.6 Example 1.

**R implementation:**

- `add_interviews()` in `R/creel-design.R`: computes
  `.expansion = n_counted / n_interviewed` for each interview row. When
  `n_counted = 0` and `n_interviewed = 0` (zero-effort site), `.expansion` is set to
  `NA`; the effort estimator treats this as zero contribution.

---

## 3. Effort Estimator (Jones & Pollock Eq. 19.4)

**Published source:** Jones & Pollock (2012), Eq. 19.4, p. 911

**Formula:**

$$\hat{E} = \sum_{i=1}^{n} \frac{e_i}{\pi_i}$$

where $e_i$ is the enumeration-expanded effort for interview $i$ (see Section 4 below),
$\pi_i$ is the inclusion probability, and the sum runs over all interview records.

**Statistical meaning:** This is a Horvitz-Thompson (HT) total estimator. Dividing by
$\pi_i$ is the HT inverse-probability weight: if a unit is sampled with probability
$\pi_i$, it represents $1/\pi_i$ units in the population. Summing the weighted
contributions gives an unbiased estimator of the population total (under the design).

**R implementation:**

- `estimate_effort_br()` in `R/creel-estimates-bus-route.R`, lines 21–215:
  - Line 69: `interviews$.e_i <- interviews[[effort_col]] * interviews$.expansion`
    (computes expanded effort; see Section 4)
  - Line 83: `interviews$.contribution <- interviews$.e_i / interviews$.pi_i`
    (computes $e_i / \pi_i$ for each row)
  - Line 96: `total_estimate <- sum(interviews$.contribution, na.rm = TRUE)`
    (sums to produce $\hat{E}$)
- `estimate_effort()` in `R/creel-estimates.R`, line 310: dispatches to
  `estimate_effort_br()` when `design$design_type == "bus_route"`.

---

## 4. Expanded Effort per Interview (eᵢ)

**Published source:** Malvestuto (1996), Box 20.6, p. 614 (implicit in the expansion step)

**Formula:**

$$e_i = \text{hours\_fished}_i \times \text{expansion}_i$$

**Statistical meaning:** The effort recorded in an interview represents only the fishing
party interviewed. Multiplying by the expansion factor scales that effort to represent
all $n\_counted$ parties at the site, not just the $n\_interviewed$ parties.

**R implementation:**

- `estimate_effort_br()` in `R/creel-estimates-bus-route.R`:
  - Line 69: `interviews$.e_i <- interviews[[effort_col]] * interviews$.expansion`

When expansion = 1 (all parties interviewed, as in Box 20.6 Example 1), $e_i$ equals the
raw effort. For example, Site C has 6 interviews each with `hours_fished = 57.5 / 6`:

$$e_C = (57.5/6) \times (6/6) = 57.5 \text{ h per interview row}$$

Wait — more precisely, each of the 6 Site C rows contributes $e_i = (57.5/6) \times 1$,
and these sum to 57.5 h before the $1/\pi_C = 5$ weight is applied, yielding the site
contribution of 287.5 angler-hours.

---

## 5. Harvest Estimator (Jones & Pollock Eq. 19.5)

**Published source:** Jones & Pollock (2012), Eq. 19.5, p. 912

**Formula:**

$$\hat{H} = \sum_{i=1}^{n} \frac{h_i}{\pi_i}$$

where $h_i = \text{harvest}_i \times \text{expansion}_i$ is the enumeration-expanded
harvest for interview $i$.

**Statistical meaning:** Structurally identical to the effort estimator (Eq. 19.4) with
harvest substituted for effort. Both are HT totals over the same sampling design.

**R implementation:**

- `estimate_harvest_br()` in `R/creel-estimates-bus-route.R`:
  computes `.h_i = harvest * .expansion` and `.contribution = .h_i / .pi_i`.
- `estimate_total_harvest()` sums those contributions via
  `survey::svytotal(~.contribution, ...)` to give $\hat{H}$.

---

## 5a. Harvest Rate (HPUE) as a Ratio of HT Totals

Eq. 19.4 and Eq. 19.5 are **totals**. Jones & Pollock define no rate estimator for
the bus-route design, so 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 the ratio-of-means form, the same quantity `estimate_harvest_rate()`
returns for standard designs, and it carries the same `method` string
(`"ratio-of-means-hpue"`). Because $e_i$ is built from angler-effort, the result
is fish per **angler**-hour, matching the CPUE denominator on the same design.

**Variance.** Computed with `survey::svyratio()` over both contribution columns,
not by dividing two separately estimated totals. $\hat{H}$ and $\hat{E}$ come from
the same interviews and are strongly positively correlated; propagating their SEs
as if independent overstates the standard error substantially (roughly eightfold
on the package's bus-route harvest fixture).

## 5b. Incomplete trips: truncated mean of ratios

Everything above assumes completed trips. Anglers intercepted mid-trip are a
different estimation problem, and `use_trips = "incomplete"` uses a different
estimator for it.

Hoenig, Jones, Pollock, Robson & Wade (1997) analyse both candidates under the
roving-type sampling that applies here — encounter probability proportional to
time fished, catch observed only up to the interview moment. 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

$$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^{*}$. That is not the population catch rate, 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 the correct expectation — the ratio of total catch to total effort.

Two departures from the paper's plain average are forced by the bus-route
design. Interviews are not equally likely, so each is weighted by
$w_i = \text{expansion}_i / \pi_i$, giving 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}$$

where $a_i$ is angler-effort, so the result is fish per **angler**-hour rather
than fish per party-hour. It is reported as `method = "mean-of-ratios-hpue"`.

**Truncation is not optional.** The mean-of-ratios estimator has *infinite*
asymptotic variance: $E[1/L_j]$ is infinite as trip length approaches zero.
Hoenig et al. recommend discarding trips shorter than 30 minutes, which is the
`truncate_at = 0.5` default. The threshold applies to elapsed trip duration, not
to angler-hours — it is the short clock interval that makes the reciprocal
explode, and a party of five fishing twelve minutes supplies a full angler-hour
while still being the unstable case. `truncate_at = NULL` disables truncation
and warns.

**Variance.** As in 5a, `survey::svyratio()` linearises over numerator and
denominator together; here the denominator is the sum of the weights.

**Comparing the two.** `use_trips = "diagnostic"` returns both slots. They use
different estimators, because each is the estimator its trip type supports, but
both report fish per angler-hour — so the gap between them is attributable to
trip status rather than to a change of physical quantity.

**R implementation:**

- `br_harvest_rate_estimates()` in `R/creel-estimates-bus-route.R` (shared ratio
  machinery for both trip paths).
- `br_incomplete_harvest_rate()` in `R/creel-estimates-bus-route.R`.
- `estimate_harvest_rate()` in `R/creel-estimates.R`: dispatches to
  `estimate_harvest_br()` when `design$design_type == "bus_route"`.

---

## 6. Variance Estimation

**Published source:** Horvitz & Thompson (1952); implemented via Taylor linearization
in the R `survey` package (Lumley 2010).

**Method:** tidycreel uses `survey::svytotal(~.contribution, svy_br)` where
`svy_br` is an `svydesign` object constructed from the interview data. The `.contribution`
column holds $e_i / \pi_i$ (or $h_i / \pi_i$). The `survey` package applies Taylor
linearization to compute the standard error of this total.

**R implementation:**

- `estimate_effort_br()`, lines 99–112: constructs `svy_br` via
  `survey::svydesign(ids = ~1, ...)`, applies `get_variance_design()` for the
  selected variance method, and calls `survey::svytotal(~.contribution, svy_br)`.
- Phase 26 cross-validation (`tests/testthat/test-cross-validation.R`) confirms that
  the `survey::svytotal` variance matches a manually constructed `svydesign` calculation
  to tolerance 1e-6.

Bootstrap (`variance = "bootstrap"`) and jackknife (`variance = "jackknife"`) alternatives
are available via `get_variance_design()`, consistent with all other tidycreel estimators.

---

## 7. Summary Traceability Table

| Quantity | Formula | Source | Page | R Location |
|----------|---------|--------|------|------------|
| Inclusion probability | πᵢ = p_site × p_period | Jones & Pollock (2012) | p. 912 | `creel-design.R`: `creel_design()`, `add_interviews()` |
| Enumeration expansion | expansion = n_counted / n_interviewed | Malvestuto (1996) Box 20.6 | p. 614 | `creel-design.R`: `add_interviews()` |
| Expanded effort | eᵢ = angler-hours × expansion | Malvestuto (1996) Box 20.6 | p. 614 | `creel-estimates-bus-route.R`: `estimate_effort_br()` |
| HT effort total | Ê = Σ(eᵢ/πᵢ) | Jones & Pollock (2012) Eq. 19.4 | p. 911 | `creel-estimates-bus-route.R` lines 83, 96 |
| Expanded harvest | hᵢ = harvest × expansion | Malvestuto (1996) Box 20.6 | p. 614 | `creel-estimates-bus-route.R` (harvest branch) |
| HT harvest total | Ĥ = Σ(hᵢ/πᵢ) | Jones & Pollock (2012) Eq. 19.5 | p. 912 | `creel-estimates-bus-route.R`: `estimate_total_harvest_br()` |
| Harvest rate (HPUE) | Ĥ / Ê | ratio of Eq. 19.5 to Eq. 19.4 | pp. 911–912 | `creel-estimates-bus-route.R`: `br_harvest_rate_estimates()` |
| Variance | Taylor linearization on Σ(eᵢ/πᵢ) | Lumley (2010) | — | `survey::svytotal(~.contribution, svy_br)` |

---

## 8. Why πᵢ Matters: A Quantitative Example

Some implementations use a fixed value such as $\pi_i = 0.5$ for all sites, or compute
$\pi_i$ from interview timing data (e.g., wait time / circuit time). Both approaches are
statistically incorrect: neither is the inclusion probability of the sampling design.

The bias can be large and heterogeneous across sites. Using the Malvestuto (1996) Box 20.6
data, here is the effect of substituting $\pi_i = 0.5$ for the correct design-based
values:

| Site | Correct πᵢ | eᵢ (h) | Correct eᵢ/πᵢ | Incorrect (π=0.5) eᵢ/πᵢ | Error |
|------|-----------|--------|--------------|-------------------------|-------|
| A    | 0.150     | 30.0   | 200.0        | 60.0                    | −70%  |
| B    | 0.125     | 20.0   | 160.0        | 40.0                    | −75%  |
| C    | 0.200     | 57.5   | 287.5        | 115.0                   | −60%  |
| D    | 0.025     |  5.0   | 200.0        | 10.0                    | −95%  |
| **Total** | | | **847.5** | **225.0** | **−73%** |

With $\pi_i = 0.5$, the total effort estimate would be 225.0 angler-hours — a 73%
underestimate of the correct 847.5. The bias is not uniform: Site D is underestimated
by 95% because it has the lowest correct $\pi_i$ (0.025) but the incorrect formula
assigns it the same weight as Site C (π = 0.5).

The direction and magnitude of bias depend entirely on the distribution of site
probabilities in the actual design. The only way to avoid bias is to use the
design-specified inclusion probabilities — which is what tidycreel does.

For the specific Site C example cited in the design documents: correct contribution
= 57.5 / 0.20 = **287.5**; with π = 0.5 it would be 57.5 / 0.5 = **115.0** — a 2.5×
underestimate for that site alone.

---

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

- Horvitz, D. G., & Thompson, D. J. (1952). A generalization of sampling without
  replacement from a finite universe. *Journal of the American Statistical Association*,
  47(260), 663–685.

- Jones, C. M., & Pollock, K. H. (2012). Recreational survey methods: estimation of effort,
  harvest, and abundance. 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.

- Malvestuto, S. P. (1996). Sampling the recreational angler. Chapter 20 in *Fisheries
  Techniques* (2nd ed.), pp. 591–623. American Fisheries Society.

- Malvestuto, S. P., Davies, W. D., & Shelton, W. L. (1978). An evaluation of the roving
  creel survey with nonuniform probability sampling. *Transactions of the American Fisheries
  Society*, 107(2), 255–262.