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:
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.
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 |
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.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:
interviews$.e_i <- interviews[[effort_col]] * interviews$.expansion
(computes expanded effort; see Section 4)interviews$.contribution <- interviews$.e_i / interviews$.pi_i
(computes \(e_i / \pi_i\) for each
row)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".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:
interviews$.e_i <- interviews[[effort_col]] * interviews$.expansionWhen 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.
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}\).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).
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".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).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.
| 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) |
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.
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.