---
title: "GLMM Aerial Effort Estimation"
output:
  rmarkdown::html_vignette:
    highlight: null
vignette: >
  %\VignetteIndexEntry{GLMM Aerial Effort Estimation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## When to Use the GLMM Estimator

If your pilot always flies at the same time of day — say, 10 AM — your
instantaneous count systematically over- or under-represents total daily effort.
Early-morning flights catch fewer anglers than are present at peak hours;
late-afternoon flights may miss early starters entirely. When flight timing is
non-random, the simple expansion (count × h_open / v) inherits this temporal
bias.

The GLMM approach (Askey et al. 2018) models how angler counts change through
the day using a quadratic hour effect and a day-level random intercept:
`count ~ poly(hour, 2) + (1 | date)`. Once the diurnal curve is estimated, the
model integrates predicted counts across the full open-water window — correcting
for wherever in the curve the actual flight fell.

**When to use the GLMM estimator:**

1. Flights always occur in the same part of the day (fixed morning or afternoon
   schedule).
2. You have multiple overflights per day across several days (minimum 8–10
   survey days recommended for stable random-effect estimation).

**When to use the simple estimator:**

1. Flight timing is randomly assigned across the open-water window.
2. You have only one count per day.

For the basic aerial workflow without GLMM correction, see the
[aerial surveys vignette](aerial-surveys.html).

## Example Data

The `example_aerial_glmm_counts` dataset contains 12 survey days with 4
overflights per day at fixed hours (7, 10, 13, and 16 hours), producing
48 observations. Angler counts follow a diurnal curve with day-level Poisson
variability — representative of a scenario where a fixed morning-to-afternoon
flight schedule is used throughout the season.

```{r data}
library(tidycreel)
data(example_aerial_glmm_counts)
head(example_aerial_glmm_counts)
```

The four columns are `date`, `day_type`, `n_anglers` (instantaneous count), and
`time_of_flight` (decimal hour of each overflight).

## Building the Aerial Design

Build an aerial `creel_design` from the survey dates and attach the count data.
The `h_open = 14` argument specifies the number of hours the fishery is open
each day — this enters the final expansion after the diurnal curve is
integrated.

```{r design}
aerial_cal <- unique(example_aerial_glmm_counts[, c("date", "day_type")])
aerial_cal <- aerial_cal[order(aerial_cal$date), ]

design <- creel_design(
  aerial_cal,
  date        = date,
  strata      = day_type,
  survey_type = "aerial",
  visibility_correction = "none",
  angler_ratio = 1,
  angler_ratio_se = 0,
  h_open      = 14
)

design <- add_counts(design, example_aerial_glmm_counts, count_col = n_anglers)
print(design)
```

## GLMM Effort Estimation

Call `estimate_effort_aerial_glmm()` with `time_col = time_of_flight`. The
default model fits a negative-binomial GLMM with a quadratic temporal effect
and a day-level random intercept:

```{r glmm-default}
glmm_result <- estimate_effort_aerial_glmm(design, time_col = time_of_flight)
print(glmm_result)
```

The model fits the diurnal count curve over all 48 observations, then
numerically integrates the predicted mean count across the full open-water
window (from 0.5 hours before the earliest flight to 0.5 hours after the
latest). The integrated area is scaled by `h_open / visibility_correction` to
convert counts to angler-hours.

## Variance and Confidence Intervals

Two variance methods are available.

**Delta method (default):** Propagates the fixed-effect covariance matrix from
`lme4::vcov()` to the derived integral via a gradient vector. This is fast
and analytic, and is the default when `boot = FALSE`.

**Parametric bootstrap:** `lme4::bootMer()` re-fits the model under parametric
resampling `nsim` times and uses the SD of the resulting totals as the SE.
This method can give more accurate CIs for skewed count distributions. Use it
for final production analyses; the delta method is appropriate for exploratory
work.

**Whichever method you choose, this vignette's design reports no interval at
all.** `creel_design()` above sets `visibility_correction = "none"`, which
declares that no detection-probability study was done. That is a statement of
ignorance, not a value of one: the correction's uncertainty is unknown, so it
cannot be propagated, and `se`, `ci_lower` and `ci_upper` all come back `NA`
rather than pretending the unknown component contributes zero.

To get an interval, say what the correction's uncertainty is — pass a numeric
`visibility_correction` together with `visibility_se` to `creel_design()`.
Declaring it known and exactly certain (`visibility_se = 0`) is also a valid
claim and does produce an interval; it is simply a different claim from having
never measured it.

```{r glmm-boot, eval = FALSE}
# Bootstrap CIs — use nboot = 500 for production analyses
glmm_boot <- estimate_effort_aerial_glmm(
  design,
  time_col = time_of_flight,
  boot = TRUE,
  nboot = 100L
)
print(glmm_boot)
```

## Downstream Estimation

`estimate_effort_aerial_glmm()` returns an effort estimate; it does not write
that estimate back into the design, and the total estimators do not take one as
an argument — `estimate_total_catch()` derives effort itself, from the counts
attached to whatever design it is given. So the GLMM figure is a standalone
result to read, report, or combine by hand, not an input the rest of the
pipeline picks up automatically.

What the downstream estimators need is a design carrying interviews. The
example below builds one from the complementary aerial interview data and
estimates catch rate and total catch from it. Note that this is a *different*
design object from the one fitted above: it uses `example_aerial_counts`, and
its effort comes from the standard aerial path rather than from the GLMM.

```{r downstream}
# Build a complementary design with matching interview dates
aerial_int_cal <- unique(example_aerial_counts[, c("date", "day_type")])
aerial_int_cal <- aerial_int_cal[order(aerial_int_cal$date), ]

design_int <- creel_design(
  aerial_int_cal,
  date        = date,
  strata      = day_type,
  survey_type = "aerial",
  visibility_correction = "none",
  angler_ratio = 1,
  angler_ratio_se = 0,
  h_open      = 14
)
design_int <- add_counts(design_int, example_aerial_counts)
design_int <- add_interviews(design_int, example_aerial_interviews,
  catch       = walleye_catch,
  effort      = hours_fished,
  trip_status = trip_status
)
catch_rate <- estimate_catch_rate(design_int)
print(catch_rate)

total_catch <- estimate_total_catch(design_int)
print(total_catch)
```

`example_aerial_interviews` carries no party-size column, so `add_interviews()`
cannot normalise effort to angler-hours and `estimate_total_catch()` warns that
the CPUE it multiplies is per *party*-hour while the count-derived effort is
angler-hours. The warning is correct for this dataset and is left visible rather
than suppressed: with one angler per party the total is right as printed, and
with larger parties it is **overstated** — by the mean party size, because the
party-hour rate is multiplied by an effort that already counts every angler.
Supplying `n_anglers = 2` on this example halves the total, from 250.6 to 125.3. Supply `n_anglers` from your own
interview data, or state a constant party size (`n_anglers = 1`) when every
interview really is a single angler, to remove the ambiguity.

## Comparison: Simple vs. GLMM Estimator

The two estimators need the counts at different grains, so they take different
designs built from the same data.

The GLMM consumes each overflight individually — the flight times are what it
fits the diurnal curve against, so the four flights a day must stay four rows.

The simple estimator does not model time. It treats each count as an
instantaneous estimate of the anglers present and expands by `h_open / v`. Four
flights on one day are therefore four looks at that day, not four sampled days:
they average to the day's mean occupancy, and the spread between them becomes
the within-day variance component. Naming `count_time_col` is what performs that
aggregation. Without it the four counts would be summed as separate sampling
units and the estimate would come back four times too large.

```{r comparison}
# Daily means for the simple estimator: the flights are four looks at each day,
# so they aggregate rather than accumulate.
design_daily <- add_counts(
  creel_design(
    aerial_cal,
    date        = date,
    strata      = day_type,
    survey_type = "aerial",
    visibility_correction = "none",
    angler_ratio = 1,
    angler_ratio_se = 0,
    h_open      = 14
  ),
  example_aerial_glmm_counts,
  count_col      = n_anglers,
  count_time_col = time_of_flight
)

# Simple aerial estimator — no diurnal correction
simple_result <- estimate_effort(design_daily)

# GLMM result from above
# glmm_result already computed

# Side-by-side comparison. Both estimators report a total across the sampled
# days, so the two numbers answer the same question and the gap between them is
# the diurnal correction.
comparison <- data.frame(
  method = c("GLMM", "Simple"),
  estimate = c(
    glmm_result$estimates$estimate,
    simple_result$estimates$estimate
  ),
  target = c(
    glmm_result$effort_target,
    simple_result$effort_target
  ),
  stringsAsFactors = FALSE
)

print(comparison)
```

The `target` column is worth checking rather than assuming. Both estimators
report `sampled_days`, so the two totals cover the same twelve days and the gap
between them is the diurnal correction and nothing else. The GLMM corrects for
the fact that all flights occurred at fixed hours (7, 10, 13, 16); the simple
estimator treats each count as representative of the full open-water window,
which inflates or deflates the total depending on where the peak falls in the
diurnal curve. Here it inflates, so the corrected estimate is the lower of the
two.

Neither estimator reports a confidence interval here, for the reason given under
"Variance and Confidence Intervals" above: this design declares
`visibility_correction = "none"`.

## Custom Formula

For surveys where a linear temporal term is sufficient — or where the analyst
prefers to control the model structure directly — pass a custom `formula`. The
formula must reference the actual count column (`n_anglers`) and the time
column by its exact name in the data.

```{r custom-formula}
glmm_linear <- estimate_effort_aerial_glmm(
  design,
  time_col = time_of_flight,
  formula  = n_anglers ~ time_of_flight + (1 | date)
)
print(glmm_linear)
```

A linear temporal term reduces flexibility but can improve stability when only
a few survey days are available. Use the default quadratic formula when you
have 8 or more survey days.

## References

- Askey, P. J., Ward, H., Godin, T., Boucher, M., & Northrup, S. (2018).
  Angler effort estimates from instantaneous aerial counts: use of
  high-frequency time-lapse camera data to inform model-based estimators.
  *North American Journal of Fisheries Management*, 38(1), 194–209.
  [https://doi.org/10.1002/nafm.10010](https://doi.org/10.1002/nafm.10010)

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