---
title: "Choosing a classifier threshold"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Choosing a classifier threshold}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
# The four fits take about half an hour on 16 cores, so the results were
# saved by dev/precompute-threshold.R and are read from
# threshold-results.rds; the refits with fp ~ temp + rh, by
# dev/precompute-threshold-fp.R, from threshold-fp-results.rds.
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 5, out.width = "100%")
library(crywolf)
sweep <- readRDS("threshold-results.rds")
sweep_fp <- readRDS("threshold-fp-results.rds")
val <- function(par, th, digits = 2, x = sweep) {
  r <- x[x$parameter == par & x$threshold == th, ]
  sprintf(paste0("%.", digits, "f (%.", digits, "f to %.", digits, "f)"),
          r$mean, r$lower, r$upper)
}
pct <- function(par, th) {
  r <- sweep[sweep$parameter == par & sweep$threshold == th, ]
  sprintf("%.2f%% (%.2f%% to %.2f%%)", 100 * r$mean, 100 * r$lower, 100 * r$upper)
}
rate <- function(th) {
  sprintf("%.1f%%", 100 * unique(sweep$detection_rate[sweep$threshold == th]))
}
```

```{r stale, echo = FALSE, results = "asis"}
cat(crywolf:::.stale_note(sweep, "threshold-results.rds"))
cat(crywolf:::.stale_note(sweep_fp, "threshold-fp-results.rds"))
```

A classifier gives each recording a confidence score, and the analyst picks
the score above which a recording counts as a detection. A low threshold
keeps more true detections and more false ones; a high threshold loses
both. The usual way to choose it is to measure precision and recall on a
validation set and take, for instance, the threshold that maximizes their
harmonic mean, F1 (Knight et al., 2017). Piatti et al. (2026) showed on the
data used here that the choice can change the estimated environmental
effects on occupancy.

A false-positive occupancy model changes the question. It estimates the
false-positive rate instead of assuming it is zero, so a threshold that
lets in more false positives should, in principle, be absorbed by a higher
estimated false-positive rate rather than by biased occupancy effects.
`fp_threshold_sweep()` checks whether that holds: it builds the occasions at
several thresholds with `fp_aggregate()`, fits the same model to each, and
collects the estimates.

## Thresholds from reviewed clips

Numeric thresholds remain available. Alternatively, pass reviewed clips from
`fp_read_validation()` to select MaxF1 and MaxPrecision automatically:

```{r automatic-thresholds, eval = FALSE}
reviews <- fp_read_validation("validation_by_clip.csv")
selected <- fp_threshold_sweep(
  occasions, detections, sites = nights, site = "site_night",
  thresholds = c("MaxF1", "MaxPrecision"), validation = reviews,
  random_recorder = FALSE
)
attr(selected, "threshold_selection")
attr(selected, "validation_metrics")
```

Each distinct reviewed score is a candidate threshold; scores equal to the
threshold count as detections. Precision measures the fraction of retained
clips that are true positives. Recall measures the fraction of reviewed
positives retained, and F1 combines precision and recall. MaxPrecision tends
to pick a high score with few retained clips, which makes it unstable in
small validation samples; look at `validation_metrics` before trusting it.
Ties favour greater
recall, then the smaller threshold. If both methods select the same threshold,
the model is fitted once.

The validation sample should represent the target species and the score range
of interest. Recall here refers to positives in that sample: it does not count
vocalizations missed by the classifier or excluded by a previous score cutoff.
Missing scores or verdicts are dropped with a warning; at least one reviewed
positive is needed. Use `species` to select the target when the validation
contains several species. Validation is used only to select thresholds; it is
not automatically used as confirmation data in the occupancy model.

## Running the sweep

The data are the one-minute recordings of the `curiango` dataset, scored by
the custom classifier of Piatti et al. (2026). Each recording is an
occasion; nights at a site are the sites of the model.

```{r sweep, eval = FALSE}
zs <- function(v) as.numeric(scale(v))
m <- curiango$minutes
nt <- curiango$nights
occasions <- data.frame(file = paste(m$site_night, format(m$time, "%Y%m%d_%H%M")),
                        site_night = m$site_night, time = m$time,
                        bio = zs(m$bio), aci = zs(m$aci), recorder = m$site)
detections <- data.frame(file = occasions$file, confidence = m$conf_custom)
sites <- data.frame(site_night = nt$site_night, temp = zs(nt$temp),
                    rh = zs(nt$rh), moon = zs(log(nt$moon)))
verified <- data.frame(site_night = curiango$verified$site_night,
                       y = curiango$verified$y)

sweep <- fp_threshold_sweep(occasions, detections,
                            thresholds = c(0.2, 0.4, 0.6, 0.8),
                            sites = sites, verified = verified,
                            site = "site_night", recorder = "recorder",
                            occ = ~ temp + rh + moon, det = ~ bio + aci, fp = ~ 1,
                            threads_per_chain = 4)
```

The function does not choose the thresholds. Here they span the range of
the scores; in a real study, include the threshold chosen on the
validation set (0.40 was the F1-maximizing threshold of this classifier in
Piatti et al., 2026) and a few on either side. The four fits took
`r round(attr(sweep, "minutes"))` minutes in all.

## Results

```{r detection-rate}
unique(sweep[c("threshold", "detection_rate")])
```

From the lowest to the highest threshold, the share of one-minute
recordings counted as detections falls from `r rate(0.2)` to `r rate(0.8)`.

```{r plot, fig.cap = "Posterior means and 95% intervals of the occupancy effects, the detection intercept, and the mean false-positive rate at four thresholds of the custom classifier."}
plot(sweep, parameters = c("occ_temp", "occ_rh", "occ_moon", "p10_bar",
                           "det_(Intercept)"))
```

```{r table}
keep <- c("occ_temp", "occ_rh", "occ_moon", "p10_bar")
tab <- sweep[sweep$parameter %in% keep,
             c("threshold", "parameter", "mean", "lower", "upper", "max_rhat")]
tab[c("mean", "lower", "upper")] <- signif(tab[c("mean", "lower", "upper")], 3)
tab$max_rhat <- round(tab$max_rhat, 3)
tab[order(tab$parameter, tab$threshold), ]
```

The mean false-positive rate per recording falls from
`r pct("p10_bar", 0.2)` at a threshold of 0.2 to `r pct("p10_bar", 0.8)` at
0.8: the model attributes many of the extra detections at low thresholds to
false positives, as intended.

The occupancy effects, however, are not stable over the whole range. From
0.4 to 0.8 they stay close: relative humidity goes from
`r val("occ_rh", 0.4)` to `r val("occ_rh", 0.8)`, and temperature from
`r val("occ_temp", 0.4)` to `r val("occ_temp", 0.8)`. At 0.2 they change:
relative humidity is `r val("occ_rh", 0.2)`, and the temperature effect
changes sign, to `r val("occ_temp", 0.2)`. A constant false-positive rate
does not absorb what a low threshold lets in.

```{r divergences}
unique(sweep[c("threshold", "divergences", "max_rhat")])
```

## False positives that follow the conditions of the night

One reading of the change at 0.2 is that the false positives let in by a
low threshold are themselves related to the conditions of the night. If so,
a constant false-positive rate assigns part of them to occupancy, and the
occupancy effects of those conditions are biased. The remedy suggested in
`vignette("crywolf")` is to model the false-positive rate with the
occupancy covariates. The thresholds 0.2 and 0.4 were refitted with
`fp = ~ temp + rh` (`dev/precompute-threshold-fp.R`; the two fits took
`r round(attr(sweep_fp, "minutes"))` minutes):

```{r fp-covariates}
pars <- c("occ_temp", "occ_rh", "fp_temp", "fp_rh", "p10_bar")
both <- rbind(
  data.frame(fp = "~ 1", sweep[sweep$threshold %in% c(0.2, 0.4) &
                                 sweep$parameter %in% pars, ]),
  data.frame(fp = "~ temp + rh", sweep_fp[sweep_fp$parameter %in% pars, ]))
both[c("mean", "lower", "upper")] <- signif(both[c("mean", "lower", "upper")], 3)
both <- both[order(both$threshold, both$parameter, both$fp),
             c("threshold", "parameter", "fp", "mean", "lower", "upper")]
rownames(both) <- NULL
both
```

At 0.2 the false-positive rate rises with temperature,
`r val("fp_temp", 0.2, x = sweep_fp)` on the logit scale, and falls slightly
with humidity, `r val("fp_rh", 0.2, x = sweep_fp)`. With these terms, the
occupancy effect of temperature shrinks from `r val("occ_temp", 0.2)` to
`r val("occ_temp", 0.2, x = sweep_fp)`, and its interval now just includes
zero. It does not return to the negative value found at 0.4 and above. The
humidity effect hardly changes, from `r val("occ_rh", 0.2)` to
`r val("occ_rh", 0.2, x = sweep_fp)`, and stays weaker than at 0.4. At 0.4
the false-positive rate falls with humidity,
`r val("fp_rh", 0.4, x = sweep_fp)`, yet the occupancy effects barely move:
humidity goes from `r val("occ_rh", 0.4)` to
`r val("occ_rh", 0.4, x = sweep_fp)`, and temperature stays at
`r val("occ_temp", 0.4, x = sweep_fp)`. Both refits had no
divergent transitions and a largest R-hat of
`r sprintf("%.3f", max(sweep_fp$max_rhat))`.

So false positives that follow the temperature of the night explain part of
the change at the low threshold, but not all of it: what remains may come
from false positives related to conditions not in the formula, or from
changes in detection, which also depends on the threshold. This is
consistent with the simulations behind crywolf (Provete et al., in preparation), in which
putting the occupancy covariates in the false-positive formula reduced the
bias of the occupancy effects without removing it. The estimates at 0.2
should therefore still not be trusted, and the stability from 0.4 upwards,
which includes the F1-maximizing threshold, is the reassuring part of the
sweep.

## Reading a sweep

* **Stable occupancy effects across thresholds** suggest that the
  false-positive part of the model is doing its job: the threshold mostly
  changes the estimated detection and false-positive rates.
* **Occupancy effects that drift with the threshold** mean that the
  threshold affects detections in a way the model does not capture, for
  instance false positives that depend on a covariate. Try a false-positive
  formula with the occupancy covariates (`vignette("crywolf")`) before
  trusting any single threshold.
* **Very high thresholds** leave few detections, and the estimates become
  wide and lean on the priors. Very low ones make the false-positive rate
  large, and the gap between true and false detection small.
* Check the diagnostics of each fit (the `divergences` and `max_rhat`
  columns) before comparing them.

## References

Knight, E. C., Hannah, K. C., Foley, G. J., Scott, C. D., Brigham, R. M., &
Bayne, E. (2017). Recommendations for acoustic recognizer performance
assessment with application to five common automated signal recognition
programs. *Avian Conservation and Ecology*, 12, 14.
https://doi.org/10.5751/ACE-01114-120214

Piatti, L., de Sousa, D. L. H., Oliveira, B. S., de Melo, A. V., Sugai,
L. S. M., & Provete, D. B. (2026). Methodological choices influence
ecological inference in passive acoustic monitoring of a Neotropical
nightjar. *Journal of Avian Biology*, in press.

Provete, D. B., et al. (in preparation).