creel survey, fisheries, design-based inference, R, tidycreel, Quarto
Chapters 1–3 established the why — why creel surveys exist, what they estimate, and what makes the estimators statistically defensible. This chapter establishes the how before the theory goes deeper.
A complete tidycreel analysis has four steps: load the data, build a design object that encodes the survey structure, run the estimators, and summarize the results. This chapter walks through all four on the Harlan Reservoir dataset and produces defensible estimates of effort, harvest, and catch rate by the end. No theory not already in Chapters 1–3; just the pipeline running on real data.
If you want to see the complete sequence in one view, it fits in about forty lines:
# ── 1. Load ──────────────────────────────────────────────────────────────────
schedule <- readRDS(here("data/harlan_schedule.rds"))
counts <- readRDS(here("data/harlan_counts.rds"))
interviews <- readRDS(here("data/harlan_interviews.rds"))
catch <- readRDS(here("data/harlan_catch.rds"))
# ── 2. Pre-process ───────────────────────────────────────────────────────────
interviews <- interviews |>
left_join(
catch |> group_by(interview_id) |>
summarise(
total_catch = sum(n_fish),
harvested = sum(n_fish[catch_type == "harvested"]),
.groups = "drop"
),
by = "interview_id"
) |>
mutate(across(c(total_catch, harvested), \(x) if_else(is.na(x), 0L, as.integer(x))))
# ── 3. Build design object ───────────────────────────────────────────────────
design <- creel_design(schedule, date = "date", strata = "day_type") |>
add_counts(counts, count_time_col = "count_time") |>
add_interviews(interviews,
catch = "total_catch", harvest = "harvested", effort = "hours_fished",
trip_duration = "hours_fished", trip_status = "trip_status",
n_anglers = "n_anglers", angler_type = "angler_type",
angler_method = "angler_method", interview_type = "access"
) |>
add_catch(catch,
catch_uid = "interview_id", interview_uid = "interview_id",
species = "species", count = "n_fish", catch_type = "catch_type"
)
# ── 4. Estimate ──────────────────────────────────────────────────────────────
effort <- estimate_effort(design, target = "stratum_total")
harvest <- estimate_total_harvest(design, target = "stratum_total")
cpue <- estimate_catch_rate(design)The rest of this chapter works through each block in detail, explaining what each function expects and what its output means. Chapters 9–16 go deeper into every step.
Four tables form the Harlan Reservoir dataset. They are linked by date and interview_id and must stay separate until the estimation step — collapsing them prematurely discards the information needed to apply the correct expansion weights. Chapter 9 covers the linked-table structure in detail; for now, looking at each table’s shape is enough.
harlan_schedule <- readRDS(here("data/harlan_schedule.rds"))
harlan_counts <- readRDS(here("data/harlan_counts.rds"))
harlan_interviews <- readRDS(here("data/harlan_interviews.rds"))
harlan_catch <- readRDS(here("data/harlan_catch.rds"))head(harlan_schedule, 3)# A tibble: 3 × 4
date reservoir day_type sampled
<date> <chr> <chr> <lgl>
1 2022-04-02 Harlan Reservoir weekend TRUE
2 2022-04-03 Harlan Reservoir weekend TRUE
3 2022-04-04 Harlan Reservoir weekday FALSE head(harlan_counts, 3)# A tibble: 3 × 9
reservoir date day_type period section count_time bank_anglers
<chr> <date> <chr> <dbl> <chr> <chr> <int>
1 Harlan Reservoir 2022-09-02 weekday 1 0 09:00:00:000 0
2 Harlan Reservoir 2022-09-02 weekday 1 0 11:15:00:000 0
3 Harlan Reservoir 2022-09-05 weekday 1 0 08:00:00:000 0
# ℹ 2 more variables: angler_boats <int>, boat_anglers <int>head(harlan_interviews, 3)# A tibble: 3 × 12
reservoir date day_type period section interview_id party_id n_anglers
<chr> <date> <chr> <dbl> <chr> <int> <chr> <dbl>
1 Harlan Res… 2022-09-01 weekday 2 0 1 6efcab0… 4
2 Harlan Res… 2022-09-01 weekday 2 0 2 d19b214… 3
3 Harlan Res… 2022-09-01 weekday 2 0 3 0093195… 3
# ℹ 4 more variables: angler_type <chr>, angler_method <chr>,
# hours_fished <dbl>, trip_status <chr>head(harlan_catch, 3)# A tibble: 3 × 4
interview_id species catch_type n_fish
<int> <chr> <chr> <int>
1 1 White Bass released 11
2 1 White Bass harvested 32
3 1 Yellow Perch released 1The four tables represent four levels of resolution: which days were in the survey and which were sampled (schedule), instantaneous counts of anglers on sampled days (counts), one record per angler party per sampled day (interviews), and one record per species per catch type per party (catch). Every downstream estimate traces back to this structure.
Before building the design, confirm the tables are internally consistent:
validate_creel_data(
counts = harlan_counts,
interviews = harlan_interviews
)── Creel Data Validation ───────────────────────────────────────────────────────63 pass | 0 warn | 0 fail── Table: counts ──✔ reservoir✔ type: class: character✔ na_rate: 0 / 236 NA (0%)✔ empty_strings: none✔ date✔ type: class: Date✔ na_rate: 0 / 236 NA (0%)✔ date_range: all within 1970-01-01 - 2100-12-31✔ day_type✔ type: class: character✔ na_rate: 0 / 236 NA (0%)✔ empty_strings: none✔ period✔ type: class: numeric✔ na_rate: 0 / 236 NA (0%)✔ negative_values: none✔ section✔ type: class: character✔ na_rate: 0 / 236 NA (0%)✔ empty_strings: none✔ count_time✔ type: class: character✔ na_rate: 0 / 236 NA (0%)✔ empty_strings: none✔ bank_anglers✔ type: class: integer✔ na_rate: 0 / 236 NA (0%)✔ negative_values: none✔ angler_boats✔ type: class: integer✔ na_rate: 0 / 236 NA (0%)✔ negative_values: none✔ boat_anglers✔ type: class: integer✔ na_rate: 0 / 236 NA (0%)✔ negative_values: none── Table: interviews ──✔ reservoir✔ type: class: character✔ na_rate: 0 / 600 NA (0%)✔ empty_strings: none✔ date✔ type: class: Date✔ na_rate: 0 / 600 NA (0%)✔ date_range: all within 1970-01-01 - 2100-12-31✔ day_type✔ type: class: character✔ na_rate: 0 / 600 NA (0%)✔ empty_strings: none✔ period✔ type: class: numeric✔ na_rate: 0 / 600 NA (0%)✔ negative_values: none✔ section✔ type: class: character✔ na_rate: 0 / 600 NA (0%)✔ empty_strings: none✔ interview_id✔ type: class: integer✔ na_rate: 0 / 600 NA (0%)✔ negative_values: none✔ party_id✔ type: class: character✔ na_rate: 0 / 600 NA (0%)✔ empty_strings: none✔ n_anglers✔ type: class: numeric✔ na_rate: 0 / 600 NA (0%)✔ negative_values: none✔ angler_type✔ type: class: character✔ na_rate: 0 / 600 NA (0%)✔ empty_strings: none✔ angler_method✔ type: class: character✔ na_rate: 0 / 600 NA (0%)✔ empty_strings: none✔ hours_fished✔ type: class: numeric✔ na_rate: 0 / 600 NA (0%)✔ negative_values: none✔ trip_status✔ type: class: character✔ na_rate: 0 / 600 NA (0%)✔ empty_strings: noneA clean pass here means no column-type errors, no impossible values, and no missing identifiers. If any check fails, validate_creel_data() returns the specific failing rows — investigate and fix before proceeding.
The creel_design object is the central data structure in tidycreel. Every estimation function takes a design object as its first argument. The object encodes the survey structure — which days were eligible, which were sampled, and at what rate — so that estimators can apply the correct expansion weights automatically. The analyst does not supply expansion factors manually; they are derived from the design.
Construction proceeds in four steps, each adding one layer of information.
creel_design() initializes the object from the sampling calendar. The date argument names the date column; strata names the column that defines sampling strata.
harlan_design <- creel_design(
calendar = harlan_schedule,
date = "date",
strata = "day_type"
)
harlan_design
── Creel Survey Design ─────────────────────────────────────────────────────────
Type: "instantaneous"
Date column: date
Strata: day_type
Calendar: 211 days (2022-04-02 to 2022-10-29)
day_type: 2 levels
Counts: "none"
Interviews: "none"
Sections: "none"At this stage the design knows the stratum structure and how many days were sampled in each stratum. The inclusion probability for each stratum — the ratio of sampled to eligible days — is implicit in these counts. Weekend days were sampled at roughly the same rate as weekdays, so their expansion weights are similar.
add_counts() attaches the count observations and identifies the column recording the time of each count. The count time matters because the instantaneous-count estimator weights each observation by its position in the sampling shift; a count at 8:00 AM represents a different portion of the day than one at 14:30.
harlan_design <- add_counts(
design = harlan_design,
counts = harlan_counts,
count_time_col = "count_time"
)Warning in svydesign.default(ids = psu_formula, strata = strata_formula, : No
weights or probabilities supplied, assuming equal probabilityharlan_design
── Creel Survey Design ─────────────────────────────────────────────────────────
Type: "instantaneous"
Date column: date
Strata: day_type
Calendar: 211 days (2022-04-02 to 2022-10-29)
day_type: 2 levels
Counts: 118 observations
PSU column: date
Count time column: count_time
Count type: "instantaneous"
Survey: <survey.design2> (constructed)
Interviews: "none"
Sections: "none"The design now holds 236 count observations linked to 118 sampled days. A survey.design2 object is constructed internally — this is the weighted survey object that estimation functions operate on.
add_interviews() requires a single column recording total catch (or harvest) per party. The Harlan interviews table does not include a precomputed catch total, so the needed column is derived from the catch table first. The join must start from harlan_interviews to preserve parties that caught nothing — they must appear as zeros, not as missing rows.
# Both columns in one pass: total_catch (all fish) and harvested (retained only)
catch_by_party <- harlan_catch |>
group_by(interview_id) |>
summarise(
total_catch = sum(n_fish),
harvested = sum(n_fish[catch_type == "harvested"]),
.groups = "drop"
)
harlan_interviews_prep <- harlan_interviews |>
left_join(catch_by_party, by = "interview_id") |>
mutate(across(c(total_catch, harvested), \(x) if_else(is.na(x), 0L, as.integer(x))))
# Confirm no rows were lost
nrow(harlan_interviews_prep) == nrow(harlan_interviews)[1] TRUETwo columns are needed: total_catch (all fish encountered, harvested + released) and harvested (fish retained). The catch_rate estimator uses total_catch; estimate_total_harvest() uses harvested. Parties with zero catch appear in harlan_interviews but not in harlan_catch; starting the join from the interviews table ensures those parties contribute zeros, not missing rows.
harlan_design <- add_interviews(
design = harlan_design,
interviews = harlan_interviews_prep,
catch = "total_catch",
harvest = "harvested",
effort = "hours_fished",
trip_duration = "hours_fished",
trip_status = "trip_status",
n_anglers = "n_anglers",
angler_type = "angler_type",
angler_method = "angler_method",
interview_type = "access"
)Warning: 179 interviews have zero catch.
ℹ Zero catch may be valid (skunked) or indicate missing data.harlan_design
── Creel Survey Design ─────────────────────────────────────────────────────────
Type: "instantaneous"
Date column: date
Strata: day_type
Calendar: 211 days (2022-04-02 to 2022-10-29)
day_type: 2 levels
Counts: 118 observations
PSU column: date
Count time column: count_time
Count type: "instantaneous"
Survey: <survey.design2> (constructed)
Interviews: 600 observations
Type: "access"
Catch: total_catch
Effort: hours_fished
Harvest: harvested
Trip status: 390 complete, 210 incomplete
Angler type: angler_type
Angler method: angler_method
Party size: n_anglers
Survey: <survey.design2> (constructed)
Sections: "none"The interview_type = "access" argument signals that these are access-point interviews — parties intercepted at a fixed exit location at the end of their trip, or occasionally mid-trip. This affects how incomplete trips are handled during estimation (Chapter 11).
add_catch() links the species-level records to the design. These are used by estimate_total_catch(), estimate_total_harvest(), and related functions that need to disaggregate catch by species.
harlan_design <- add_catch(
design = harlan_design,
data = harlan_catch,
catch_uid = "interview_id",
interview_uid = "interview_id",
species = "species",
count = "n_fish",
catch_type = "catch_type"
)
harlan_design
── Creel Survey Design ─────────────────────────────────────────────────────────
Type: "instantaneous"
Date column: date
Strata: day_type
Calendar: 211 days (2022-04-02 to 2022-10-29)
day_type: 2 levels
Counts: 118 observations
PSU column: date
Count time column: count_time
Count type: "instantaneous"
Survey: <survey.design2> (constructed)
Interviews: 600 observations
Type: "access"
Catch: total_catch
Effort: hours_fished
Harvest: harvested
Trip status: 390 complete, 210 incomplete
Angler type: angler_type
Angler method: angler_method
Party size: n_anglers
Survey: <survey.design2> (constructed)
Catch Data: 987 rows, 11 species
harvested: 430, released: 557
Sections: "none"The full design object is complete. It holds the calendar, count observations, interview records, and species-level catch — everything the estimators need.
In practice the four-step build is written as a single pipe chain. The intermediate objects are never needed separately, so there is no reason to name them:
harlan_design <- creel_design(
calendar = harlan_schedule,
date = "date",
strata = "day_type"
) |>
add_counts(
counts = harlan_counts,
count_time_col = "count_time"
) |>
add_interviews(
interviews = harlan_interviews_prep,
catch = "total_catch",
harvest = "harvested",
effort = "hours_fished",
trip_duration = "hours_fished",
trip_status = "trip_status",
n_anglers = "n_anglers",
angler_type = "angler_type",
angler_method = "angler_method",
interview_type = "access"
) |>
add_catch(
data = harlan_catch,
catch_uid = "interview_id",
interview_uid = "interview_id",
species = "species",
count = "n_fish",
catch_type = "catch_type"
)Warning in svydesign.default(ids = psu_formula, strata = strata_formula, : No
weights or probabilities supplied, assuming equal probabilityWarning: 179 interviews have zero catch.
ℹ Zero catch may be valid (skunked) or indicate missing data.harlan_design
── Creel Survey Design ─────────────────────────────────────────────────────────
Type: "instantaneous"
Date column: date
Strata: day_type
Calendar: 211 days (2022-04-02 to 2022-10-29)
day_type: 2 levels
Counts: 118 observations
PSU column: date
Count time column: count_time
Count type: "instantaneous"
Survey: <survey.design2> (constructed)
Interviews: 600 observations
Type: "access"
Catch: total_catch
Effort: hours_fished
Harvest: harvested
Trip status: 390 complete, 210 incomplete
Angler type: angler_type
Angler method: angler_method
Party size: n_anglers
Survey: <survey.design2> (constructed)
Catch Data: 987 rows, 11 species
harvested: 430, released: 557
Sections: "none"harlan_design now holds the full survey structure and is ready for estimation.
Three estimation functions cover the most common creel survey outputs: total effort, total harvest, and catch rate. All take the design object as their first argument.
estimate_effort() applies the stratified Horvitz-Thompson estimator to the count data, expanding sampled-day observations to season totals. The target = "stratum_total" argument requests the full-season expansion (as opposed to a summary of sampled days only).
effort <- estimate_effort(harlan_design, target = "stratum_total")Warning: Count variable bank_anglers contains 64 zero values.
ℹ Zero values may indicate days with no fishing activity or data collection
issues.
ℹ Consider whether zeros are true zeros or missing data.Warning: Count variable angler_boats contains 30 zero values.
ℹ Zero values may indicate days with no fishing activity or data collection
issues.
ℹ Consider whether zeros are true zeros or missing data.Warning: Count variable boat_anglers contains 30 zero values.
ℹ Zero values may indicate days with no fishing activity or data collection
issues.
ℹ Consider whether zeros are true zeros or missing data.effort$estimates# A tibble: 1 × 7
estimate se se_between se_within ci_lower ci_upper n
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 315. 9.79 9.79 0 295. 334. 118Total seasonal effort at Harlan is approximately 315 angler-hours (95% CI: 295–334). The standard error reflects between-day variability in angler counts — on some sampled days the reservoir was nearly empty, on others it was busy. That day-to-day spread, scaled by the expansion weight for each stratum, is the primary source of uncertainty in the seasonal total.
Effort broken down by stratum shows where the season’s angling time was concentrated:
effort_strata <- estimate_effort(
harlan_design,
by = day_type,
target = "stratum_total"
)Warning: Count variable bank_anglers contains 64 zero values.
ℹ Zero values may indicate days with no fishing activity or data collection
issues.
ℹ Consider whether zeros are true zeros or missing data.Warning: Count variable angler_boats contains 30 zero values.
ℹ Zero values may indicate days with no fishing activity or data collection
issues.
ℹ Consider whether zeros are true zeros or missing data.Warning: Count variable boat_anglers contains 30 zero values.
ℹ Zero values may indicate days with no fishing activity or data collection
issues.
ℹ Consider whether zeros are true zeros or missing data.effort_strata$estimates# A tibble: 2 × 8
day_type estimate se se_between se_within ci_lower ci_upper n
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 weekday 226. 8.28 8.28 0 210. 242. 83
2 weekend 88.9 5.21 5.21 0 78.6 99.2 35Figure 8.1 shows the stratum-level estimates with 95% confidence intervals.
effort_strata$estimates |>
mutate(label = c(weekday = "Weekday", weekend = "Weekend")[day_type]) |>
ggplot(aes(x = label, y = estimate, fill = day_type)) +
geom_col(alpha = 0.85) +
geom_errorbar(
aes(ymin = ci_lower, ymax = ci_upper),
width = 0.25,
linewidth = 0.7,
colour = "grey30"
) +
scale_fill_manual(
values = c(weekday = "#2171b5", weekend = "#cb181d"),
guide = "none"
) +
scale_y_continuous(labels = scales::comma) +
labs(
x = NULL,
y = "Estimated angler-hours (seasonal total)"
) +
theme_creel()
estimate_total_harvest() multiplies interview-based harvest rates by the design-weighted effort estimate to produce a season total with variance. It uses the same design object and expansion structure as estimate_effort().
harvest <- estimate_total_harvest(harlan_design, target = "stratum_total")
harvest$estimates# A tibble: 1 × 5
estimate se ci_lower ci_upper n
<dbl> <dbl> <dbl> <dbl> <int>
1 199. 16.7 167. 232. 600Estimated total harvest across the season is approximately 199 fish (95% CI: 167–232). As with effort, the confidence interval reflects the between-day variability in harvest rates across sampled days.
estimate_catch_rate() returns the ratio-of-means catch per unit effort. Unlike the totals above, CPUE does not use target = "stratum_total" because it is a rate, not a total — it is already expressed per angler-hour and does not need expansion.
cpue <- estimate_catch_rate(harlan_design)ℹ Using complete trips for CPUE estimation
(n=390, 65% of 600 interviews) [default]cpue$estimates# A tibble: 1 × 5
estimate se ci_lower ci_upper n
<dbl> <dbl> <dbl> <dbl> <int>
1 1.51 0.0928 1.32 1.69 390CPUE is 1.51 fish per angler-hour (95% CI: 1.32–1.69). This rate is computed over all species combined and all angler types. Chapter 14 disaggregates CPUE by angler type, by species, and compares the ratio-of-means estimator to the mean-of-ratios alternative.
The three estimates can be assembled into a single table. This is the minimum viable results report — a starting point for the fuller reporting format in Chapter 16.
bind_rows(
effort$estimates |>
select(estimate, se, ci_lower, ci_upper) |>
mutate(quantity = "Effort", unit = "angler-hours"),
harvest$estimates |>
select(estimate, se, ci_lower, ci_upper) |>
mutate(quantity = "Total harvest", unit = "fish"),
cpue$estimates |>
select(estimate, se, ci_lower, ci_upper) |>
mutate(quantity = "CPUE", unit = "fish / hr")
) |>
mutate(
estimate = round(estimate, 2),
se = round(se, 2),
ci_lower = round(ci_lower, 2),
ci_upper = round(ci_upper, 2)
) |>
select(quantity, unit, estimate, se, ci_lower, ci_upper)# A tibble: 3 × 6
quantity unit estimate se ci_lower ci_upper
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Effort angler-hours 315. 9.79 295. 334.
2 Total harvest fish 199. 16.7 167. 232.
3 CPUE fish / hr 1.51 0.09 1.32 1.69The CV for each estimate — se / estimate × 100 — tells you how precise the estimate is relative to its own scale. An effort CV below 20% is adequate for most management decisions; catch rate CVs are typically higher because they compound uncertainty from both the effort and interview streams.
bind_rows(
effort$estimates |> mutate(quantity = "Effort"),
harvest$estimates |> mutate(quantity = "Total harvest"),
cpue$estimates |> mutate(quantity = "CPUE")
) |>
mutate(
cv_pct = round(se / estimate * 100, 1)
) |>
select(quantity, estimate, se, cv_pct)# A tibble: 3 × 4
quantity estimate se cv_pct
<chr> <dbl> <dbl> <dbl>
1 Effort 315. 9.79 3.1
2 Total harvest 199. 16.7 8.4
3 CPUE 1.51 0.0928 6.2Harvest by species is available directly from the catch table as a raw interview summary. This is not a design-weighted estimate — it represents what the interviewed sample reported, not an expansion to the full season. Chapter 13 produces the design-weighted species-specific harvest totals.
harlan_catch |>
filter(catch_type == "harvested") |>
group_by(species) |>
summarise(total_n = sum(n_fish), .groups = "drop") |>
filter(total_n > 0) |>
mutate(species = reorder(species, total_n)) |>
ggplot(aes(x = total_n, y = species)) +
geom_col(fill = "#2171b5", alpha = 0.85) +
labs(
x = "Fish harvested (interview counts, not design-weighted)",
y = NULL
) +
theme_creel()
Figure 8.2 shows the species composition of angler harvest in the interview sample. White Bass and Walleye dominate. The species mix matters for management: a harvest regulation that targets Walleye — the most regulated species in most Nebraska reservoir fisheries — should be informed by a design-weighted Walleye-specific harvest estimate, not the raw interview count shown here.
The overall CPUE of 1.51 fish per angler-hour and total harvest of 199 fish are fishery-wide averages that pool boat and bank anglers. The interview data show whether those averages conceal important within-fishery differences.
harlan_interviews_prep |>
filter(trip_status == "complete") |>
ggplot(aes(x = angler_type, y = harvested, fill = angler_type)) +
geom_jitter(alpha = 0.25, width = 0.18, colour = "grey40", size = 0.8) +
geom_boxplot(alpha = 0.65, outlier.shape = NA, width = 0.45) +
scale_fill_manual(
values = c(bank = "#cb181d", boat = "#2171b5"),
guide = "none"
) +
scale_x_discrete(labels = c(bank = "Bank", boat = "Boat")) +
labs(
x = NULL,
y = "Fish harvested per party (complete trips)"
) +
theme_creel()
Figure 8.3 shows that boat anglers tend to harvest more fish per party than bank anglers, but both distributions are right-skewed with many zero-harvest parties. Harvest rates that differ between access types mean the pooled estimate in this chapter is a weighted average — weighted by how many angler-hours came from each type. Whether that pooled average is the right quantity for the management question depends on whether the question applies to the whole fishery or to one access type specifically. Chapter 13 demonstrates domain-restricted harvest estimation that produces separate estimates for boat and bank anglers.
This pipeline produces correct estimates for the Harlan dataset because the dataset was purpose-built to run cleanly. Real data rarely cooperates. The chapters that follow address the complications:
creel_design object stores internally and how to inspect or audit it.The pipeline here is the skeleton. The rest of the book puts flesh on it.
Four functions — creel_design(), add_counts(), add_interviews(), add_catch() — build the design object. Three functions — estimate_effort(), estimate_total_harvest(), estimate_catch_rate() — produce the estimates. Everything else in the tidycreel workflow is a variation on this structure: more estimation targets, more diagnostic checks, more domain-restricted summaries.
The numbers from this chapter — 315 angler-hours, 199 fish harvested, 1.51 fish per angler-hour — are the same numbers that appear throughout Chapters 12–16 as the estimation chapters develop each estimator in detail. Running this pipeline first means the detailed chapters are extending a result you have already seen, not building toward a result you have not.