---
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; it sets the length of the window over which 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)
```

`add_counts()` warns that the 48 rows repeat 12 sampling days with no count time
to tell them apart. For most estimators that is a mistake, and `estimate_effort()`
refuses such a design. Here it is the point: the GLMM needs each overflight as
its own row, and it models the repeats through the day-level random intercept.
So the warning is expected on a design built for `estimate_effort_aerial_glmm()`.
The comparison further down builds a second design with `count_time_col`, for
the simple estimator.

## 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 over a window `h_open` hours
long. The window starts at `open_start` if `creel_design()` was given one, and
otherwise 0.5 hours before the earliest flight. The integral is already in
people-hours, so it is scaled by the angler ratio and the visibility correction
(`angler_ratio / visibility_correction`) and not by `h_open` again.

The fit also warns `iteration limit reached`. That comes from `lme4::glmer.nb()`,
whose estimate of the negative-binomial dispersion did not settle within its
iteration limit on these 48 counts. It is a convergence warning about this
example, left visible rather than suppressed, and is open as tidycreel issue #492
along with other points about this estimator. Treat the estimate below with that in mind,
and check convergence on your own data.

## 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)
```

This example leaves out on purpose the party size that `example_aerial_interviews`
carries (`n_anglers`), to show what happens without it. `add_interviews()` warns
that it is assuming one angler per interview, and `estimate_total_catch()` warns
that the CPUE it multiplies is per *party*-hour while the count-derived effort
is angler-hours. Both warnings are correct and are left visible: with one angler
per party the total is right as printed, and with larger parties it is
**overstated**, roughly by the mean party size, because a party-hour rate is
multiplied by an effort that already counts every angler. Here the total is
3507.6 with an unknown unit; passing `n_anglers = n_anglers` gives 1838.4 fish,
as in `vignette("aerial-surveys")`. 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.

The total's standard error prints as `NA`. This design declares
`visibility_correction = "none"`: no detection study was done, so the
uncertainty of detection is unknown, and a total that leaves it out would be
a lower bound reported as an SE. The total effort's SE is `NA` for the same
reason. Supply `visibility_correction` with `visibility_se` from a detection
study to get a finite SE (see `vignette("aerial-surveys")`).

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

## Estimating Within Strata

Pass `by` to get one row per stratum, typically the design strata. The strata
enter ONE model as additive fixed effects —
`n_anglers ~ poly(time_of_flight, 2) + day_type + (1 | date)` — so weekdays and
weekends share the diurnal curve's shape and differ in its level, and each
stratum expands by its own sampled days. This follows Askey et al. (2018), where
day type is additive in every model structure compared. Fitting one model rather
than one per stratum shares strength across strata, which matters when a stratum
has few flights. The grouped path uses the delta method only (`boot = FALSE`).

The design below declares a measured visibility correction, so the grouped rows
carry standard errors; with `"none"` they would be `NA`, as explained above.

```{r glmm-by}
design_v <- creel_design(
  aerial_cal,
  date        = date,
  strata      = day_type,
  survey_type = "aerial",
  visibility_correction = 0.85,
  visibility_se = 0.05,
  angler_ratio = 1,
  angler_ratio_se = 0,
  h_open      = 14
)
design_v <- add_counts(design_v, example_aerial_glmm_counts, count_col = n_anglers)

glmm_by_day <- estimate_effort_aerial_glmm(
  design_v,
  time_col = time_of_flight,
  by = day_type
)
print(glmm_by_day)
```

The shared-shape assumption can fail: Smucker et al. (2010) report weekday and
weekend diurnal effort that differ in shape. The grouped fit therefore also fits
the time × stratum interaction and reports the BIC difference as a message. A
negative difference favours separate shapes; fit them by passing `formula`, for
example `n_anglers ~ poly(time_of_flight, 2) * day_type + (1 | date)`. The
estimate is never switched automatically, because a model chosen from the data
is not reflected in the reported standard error.

**Combining strata.** The stratum estimates are correlated: every stratum is
predicted from the same fixed effects, and the visibility correction and
angler-to-people ratio are single estimates that multiply every stratum. The
result carries their joint covariance as `strata_vcov`. Use it for the SE of the
summed total; adding the rows' SEs in quadrature treats them as independent and
understates it.

```{r glmm-by-combine}
# SE of the summed total, from the joint covariance
sqrt(sum(glmm_by_day$strata_vcov))

# The quadrature sum, shown only to compare -- it understates the SE
sqrt(sum(glmm_by_day$estimates$se^2))
```

## Results as Data

The per-stratum result converts to a tibble with its method and variance method
as columns, and plots. This design declares `visibility_correction = "none"`, so
its intervals are `NA` (see above) and the plot shows the point estimates
alone; a design with `visibility_se` draws an interval around each.

```{r results, fig.width = 6, fig.height = 3, fig.alt = "GLMM effort point estimates for weekdays and weekends; no confidence intervals, because the visibility uncertainty is unknown."}
tibble::as_tibble(glmm_by_day)

ggplot2::autoplot(glmm_by_day)
```

The `unit` column reads `NA`: the GLMM estimator does not yet record that its
effort is in angler-hours, so the plot axis carries no unit either. That is
also part of #492.

## 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 angler survey methods:
  estimation of effort, harvest, and released catch. Chapter 19 in
  *Fisheries Techniques* (3rd ed.), pp. 883–919. American Fisheries Society.

- Smucker, B. J., Lorantas, R. M., & Rosenberger, J. L. (2010). Correcting
  bias introduced by aerial counts in angler effort estimation. *North
  American Journal of Fisheries Management*, 30(4), 1051–1061.
  [https://doi.org/10.1577/M09-193.1](https://doi.org/10.1577/M09-193.1)