---
title: "tidycreel and the survey Package: A Side-by-Side Guide"
output:
  rmarkdown::html_vignette:
    highlight: null
vignette: >
  %\VignetteIndexEntry{tidycreel and the survey Package: A Side-by-Side Guide}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

For a design with instantaneous counts and access interviews, tidycreel's
effort, catch rate and total catch are built from `survey` package calls:
`svytotal()` for effort, `svyratio()` for the rate, and `svyby()` to do each
stratum separately. This page reproduces each of those three numbers with
`survey` directly, so the two can be read side by side.

It covers that one design. Other designs add steps of their own (bus-route
inclusion probabilities, aerial calibration, camera ratios), described on their
own pages; the reasoning behind each choice here is in the Explanation pages
`vignette("effort-pipeline")` and `vignette("catch-pipeline")`.

## The data

The season from `vignette("tidycreel")`: 120 days in two strata, 36 of them
sampled, one instantaneous count on each sampled day, and access interviews.

```{r data, message = FALSE}
library(tidycreel)
library(survey)

table(example_season_calendar$day_type) # days in the season, N_h
table(example_season_counts$day_type) # days sampled, n_h
head(example_season_counts, 3)
head(example_season_interviews, 3)
```

The tidycreel version, run once and used for comparison below:

```{r tidycreel, message = FALSE}
design <- creel_design(example_season_calendar, date = date, strata = day_type)
design <- add_counts(design, example_season_counts,
  count_col = total_anglers,
  period_length_col = daylight_hours
)
design <- add_interviews(design, example_season_interviews,
  catch = catch_total,
  effort = hours_fished,
  trip_status = trip_status,
  n_anglers = n_anglers
)

effort <- estimate_effort(design, target = "period_total")
rate <- estimate_catch_rate(design)
total <- estimate_total_catch(design, target = "period_total")
```

## Effort: `svytotal()`

A day's effort is its count times the length of the period the count was
randomised within: anglers present × hours open, in angler-hours. The days are
a stratified sample, and each sampled day stands for $N_h / n_h$ days of its
stratum:

```{r raw-effort}
counts <- example_season_counts
counts$angler_hours <- counts$total_anglers * counts$daylight_hours

N_h <- table(example_season_calendar$day_type)
n_h <- table(counts$day_type)
counts$w <- as.numeric(N_h[counts$day_type] / n_h[counts$day_type])

count_design <- svydesign(
  ids = ~1, strata = ~day_type, weights = ~w, data = counts
)
svytotal(~angler_hours, count_design)

effort$estimates
effort$unit
```

The estimate and standard error are the same. tidycreel's confidence interval
uses a $t$ quantile on the design's degrees of freedom, 36 days less 2 strata;
`confint()` uses the normal quantile unless it is given them:

```{r raw-effort-ci}
confint(svytotal(~angler_hours, count_design), df = degf(count_design))
```

### No finite-population correction

The design above has no finite-population correction (FPC). Survey texts would
usually add one, because days are drawn without replacement:

```{r raw-effort-fpc}
counts$fpc <- as.numeric(N_h[counts$day_type])
svytotal(~angler_hours, svydesign(
  ids = ~1, strata = ~day_type, fpc = ~fpc, data = counts
))
```

The standard error falls from about 1,286 to 1,079. tidycreel leaves the
correction out when a day has one count. The between-day variance then also
carries the variation from when in the day the count fell, and sampling more
days does not remove that part. With two or more counts on every sampled day of
a stratum, tidycreel estimates the within-day part separately and applies the
correction to the between-day part. The reasoning is in
`vignette("effort-pipeline")`; up to tidycreel 8.0.0 the correction was
always applied (#505).

## Catch rate: `svyratio()`

The rate comes from complete trips. Its denominator is **angler**-hours, each
party's hours times its number of anglers. The interviews are stratified by day
type and weighted equally:

```{r raw-rate}
trips <- subset(example_season_interviews, trip_status == "complete")
trips$angler_hours <- trips$hours_fished * trips$n_anglers

trip_design <- svydesign(
  ids = ~1, strata = ~day_type, weights = ~1, data = trips
)
svyratio(~catch_total, ~angler_hours, trip_design)

rate$estimates
rate$unit
```

This is the ratio of means (total catch over total effort), the estimator for
access interviews; `vignette("estimator-choice")` covers when tidycreel uses the
mean of ratios instead. Without `n_anglers`, the denominator would be party
hours and the rate would be per party-hour, about 0.37 here instead of 0.22.
The two are different quantities, and tidycreel reports which one it computed
in `unit`.

## Total catch: `svyby()` and a product variance

The total is not season effort times the season rate. It is built **within
each stratum** and the strata are then added, $\hat{C} = \sum_h \hat{E}_h
\hat{R}_h$, so that each stratum's rate is weighted by that stratum's effort:

```{r raw-total}
E <- svyby(~angler_hours, ~day_type, count_design, svytotal)
R <- svyby(~catch_total, ~day_type, trip_design, svyratio,
  denominator = ~angler_hours
)
E
R

e <- coef(E)
r <- coef(R)
se_e <- SE(E)
se_r <- SE(R)

sum(e * r)
sqrt(sum(e^2 * se_r^2 + r^2 * se_e^2 - se_e^2 * se_r^2))

total$estimates
total$unit
```

The variance of each stratum's product is Goodman's unbiased estimator: the
two delta-method terms minus the product of the two variances. The familiar
two-term delta method is `product_variance = "first_order"`. Both assume the
effort and rate errors are independent; `vignette("catch-pipeline")` derives the
formula and discusses that assumption.

## Mapping table

| Step | survey package | tidycreel |
|---|---|---|
| Day's effort | `count × period length` | `add_counts(count_col =, period_length_col =)` |
| Count design | `svydesign(ids = ~1, strata = ~day_type, weights = ~N_h/n_h)`, no `fpc` with one count a day | `creel_design()` + `add_counts()` |
| Effort | `svytotal(~angler_hours, ...)` | `estimate_effort(target = "period_total")` |
| Interview design | complete trips; `svydesign(ids = ~1, strata = ~day_type, weights = ~1)` | `add_interviews()` |
| Catch rate | `svyratio(~catch, ~angler_hours, ...)` | `estimate_catch_rate()` |
| Total catch | `svyby()` per stratum; $\sum_h \hat{E}_h \hat{R}_h$; Goodman variance | `estimate_total_catch(target = "period_total")` |
| Confidence interval | effort and total catch: `confint(..., df = degf(design))`; catch rate: `confint(...)`, the normal interval | `ci_lower`, `ci_upper` |

What tidycreel adds on top of these calls is the bookkeeping: the unit of each
result, the strata a total must be built within, the refusals when the data
cannot support an estimate (`vignette("refusals-and-na")`), and the bootstrap
and jackknife variance methods (`vignette("replicate-designs")`).

## The internal design object

`as_creel_svydesign()` returns the count design tidycreel built, for use with
other `survey` functions. It is the design over the **sampled days** with every
weight 1, not the season design above, and its count column now holds each
day's count × period in angler-hours, under the column's original name:

```{r internal, message = FALSE}
internal <- as_creel_svydesign(design)
head(internal$variables[, c("date", "day_type", "total_anglers")], 3)
svytotal(~total_anglers, internal)
```

That total is effort over the 36 sampled days. `estimate_effort()` expands it
to the season, and is the way to get the season total.