---
title: "Reporting a crywolf analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Reporting a crywolf analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
has_stan <- crywolf:::.has_cmdstan()
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 6.5, fig.height = 4, out.width = "100%",
                      eval = has_stan)   # model fits run only with CmdStan
# inline values come from the fit; without CmdStan they print as NA
num <- function(x, digits = 2) {
  if (has_stan) formatC(x, digits = digits, format = "f") else "NA"
}
txt <- function(x) if (has_stan) as.character(x) else "NA"
```

A Bayesian occupancy analysis is reproducible and open to scrutiny only if
the paper says what the model was, how it was computed, how well the
computation went, and how much the conclusions depend on the priors. The
Bayesian Analysis Reporting Guidelines (BARG; Kruschke, 2021) list these
items, and Davis & Kay (2023) give advice on writing statistical methods for
ecologists. This article goes through the items for a crywolf analysis and
shows where each piece of information comes from in the package. It ends
with a template for the methods and results.

We use a simulated dataset so that the article runs anywhere; with real data
every step is the same.

```{r fit, message = FALSE}
library(crywolf)
sim <- simulate_fp(nsite = 200, nocc = 10, frac_verified = 0.3, seed = 1)
m <- fp_occu(sim$data, occ = ~ x, fp = ~ 1, random_recorder = FALSE,
             threads_per_chain = 1, seed = 1, refresh = 0, show_messages = FALSE)
```

## 1. The model and why

State the design (multiple methods, Miller et al., 2011, or observation
confirmation, Chambert et al., 2015), what an occasion is, which method is
automated and which is checked, and how the checked occasions were chosen.
Then write the model, with every submodel and its covariates. `print()`
shows the formulas as fitted:

```{r print}
m
```

The equations of `vignette("crywolf")` can be copied into a paper. Say what
the gap $\Delta$ is for: it keeps the true-detection probability above the
false-positive rate at the reference values of the covariates, which rules
out the mirror-image solution in which the two are swapped.

## 2. Priors, with a reason for each

List every prior. `fp_priors()` prints the ones used:

```{r priors}
m$priors
```

Give a reason for each, even a short one: the intercept and slope priors
are weakly informative on the logit scale, and the gap prior expresses the
belief that the classifier is right more often than not. If the gap prior
came from a reviewed sample (`fp_gap_prior()`), give the counts and how the
sample was drawn. A prior predictive check (`fp_prior_predictive()`) shows
what the priors imply for data like yours, and is worth a sentence.

## 3. Computation

Report the software and versions, the number of chains, warm-up and
sampling iterations, and the target acceptance rate. These are stored in
the fit:

```{r computation}
c(crywolf = as.character(utils::packageVersion("crywolf")),
  cmdstanr = as.character(utils::packageVersion("cmdstanr")),
  CmdStan = cmdstanr::cmdstan_version())
unlist(m$sampler)
```

Then the diagnostics: the largest R-hat (Vehtari et al., 2021), the smallest
bulk and tail effective sample sizes, and the number of divergent
transitions. `summary()` gives them per parameter, and the fit keeps the
extremes:

```{r diagnostics}
unlist(m$diagnostics)
```

If the chains used more than one thread (`threads_per_chain`, `"auto"` by
default, recorded in `m$sampler`), say so: the posterior is the same, but
the draws cannot be reproduced exactly from the seed. Fit with
`threads_per_chain = 1` when they must be.

## 4. Model adequacy and comparison

A posterior predictive check tells whether the model reproduces the
features of the data that matter. Report the discrepancy measure and the
Bayesian p-value:

```{r ppc}
p <- fp_ppc(m, ndraws = 500, seed = 1)
p$p_value
```

Values near 0 or 1 flag a misfit. This model is the one that generated the
data, yet the p-value is low: a single small value calls for a look at
which sites contribute most to the discrepancy (`plot(p)`), not for
rejecting the model outright.

If models were compared, give the method (leave-one-site-out
cross-validation; Vehtari et al., 2017), the difference in expected log
predictive density with its standard error, and any sites with Pareto $k$
above 0.7, with how they were handled (`loo(reloo = TRUE)`).

## 5. The posterior

Report each quantity of interest with its posterior mean or median and a
credible interval, and say which interval (crywolf reports 95%
equal-tailed intervals). Report probabilities on the scale readers
understand, such as the false-positive rate (`p10_bar`) rather than its
logit.

```{r posterior}
s <- summary(m)$table
cols <- c("variable", "mean", "lower", "upper")
s[s$variable %in% c("occ_x", "p11_ref", "p10_bar"), cols]
```

## 6. Sensitivity to the priors

Refit with different priors and report whether the conclusions change. For
false-positive models the gap prior and the slope priors are the natural
candidates:

```{r sensitivity, message = FALSE}
refit <- function(priors) {
  f <- fp_occu(sim$data, occ = ~ x, fp = ~ 1, random_recorder = FALSE,
               priors = priors, threads_per_chain = 1, seed = 1, refresh = 0,
               show_messages = FALSE)
  t <- summary(f)$table
  t[t$variable %in% c("occ_x", "p10_bar"), cols]
}
sens <- list(
  default = s[s$variable %in% c("occ_x", "p10_bar"), cols],
  narrow_gap = refit(fp_priors(gap_mean = 2, gap_sd = 0.5)),
  tight_slopes = refit(fp_priors(slope_sd = 1))
)
out <- purrr::imap(sens, function(t, nm) data.frame(priors = nm, t)) |>
  purrr::list_rbind()
out[c("mean", "lower", "upper")] <- round(out[c("mean", "lower", "upper")], 3)
out
```

```{r max-change, include = FALSE}
slope <- out$mean[out$variable == "occ_x"]
max_change <- max(abs(slope - slope[1]))
```

Here the posterior mean of the occupancy slope moves by at most
`r num(max_change, 3)` across the three sets of priors. When they do change, report it: it means the data alone do
not settle the question.

## 7. Data and code

Share the data and the script that fits the model, with the crywolf
version. The fit object is not needed: the script and the seed reproduce it
when one thread per chain is used.

## A template

The text below is filled in from the fit above; adapt it to the study.

> We fitted a multiple-methods false-positive occupancy model (Miller et al.,
> 2011) in a Bayesian framework with the R package crywolf
> `r txt(utils::packageVersion("crywolf"))`, which uses Stan through
> CmdStan `r txt(cmdstanr::cmdstan_version())`. Occupancy probability depended on
> the standardized covariate *x*; the probability of a true automated
> detection and the false-positive rate were constant. The
> false-positive intercept was the true-detection intercept minus a positive
> gap, with a Normal(`r txt(m$priors$gap_mean)`, `r txt(m$priors$gap_sd)`) prior
> truncated at zero; intercepts had Normal(0, `r txt(m$priors$intercept_sd)`)
> priors and slopes Normal(0, `r txt(m$priors$slope_sd)`) priors on the logit
> scale. We ran `r txt(m$sampler$chains)` chains of `r txt(m$sampler$iter_warmup)`
> warm-up and `r txt(m$sampler$iter_sampling)` sampling iterations. The largest
> R-hat was `r num(m$diagnostics$max_rhat, 3)`, the smallest bulk
> effective sample size was `r txt(round(m$diagnostics$min_ess_bulk))`, and there
> were `r txt(m$diagnostics$divergences)` divergent transitions. A posterior
> predictive check based on the Freeman-Tukey discrepancy gave a Bayesian
> p-value of `r num(p$p_value)`.
>
> Occupancy `r txt(if (s$mean[s$variable == "occ_x"] < 0) "decreased" else "increased")` with *x* (slope `r num(s$mean[s$variable == "occ_x"])`,
> 95% credible interval `r num(s$lower[s$variable == "occ_x"])` to
> `r num(s$upper[s$variable == "occ_x"])`). The mean false-positive rate per
> automated occasion was `r num(s$mean[s$variable == "p10_bar"], 3)`
> (`r num(s$lower[s$variable == "p10_bar"], 3)` to
> `r num(s$upper[s$variable == "p10_bar"], 3)`). With a narrower gap prior or
> tighter slope priors, the posterior mean of the slope changed by at most
> `r num(max_change, 3)`.

## References

Chambert, T., Miller, D. A. W., & Nichols, J. D. (2015). Modeling false
positive detections in species occurrence data under different study
designs. *Ecology*, 96, 332–339. https://doi.org/10.1890/14-1507.1

Davis, A. J., & Kay, S. (2023). Writing statistical methods for ecologists.
*Ecosphere*, 14, e4539. https://doi.org/10.1002/ecs2.4539

Kruschke, J. K. (2021). Bayesian analysis reporting guidelines. *Nature
Human Behaviour*, 5, 1282–1291. https://doi.org/10.1038/s41562-021-01177-7

Miller, D. A., Nichols, J. D., McClintock, B. T., Grant, E. H. C., Bailey,
L. L., & Weir, L. A. (2011). Improving occupancy estimation when two types of
observational error occur: non-detection and species misidentification.
*Ecology*, 92, 1422–1428. https://doi.org/10.1890/10-1396.1

Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model
evaluation using leave-one-out cross-validation and WAIC. *Statistics and
Computing*, 27, 1413–1432. https://doi.org/10.1007/s11222-016-9696-4

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C.
(2021). Rank-normalization, folding, and localization: an improved R-hat
for assessing convergence of MCMC (with discussion). *Bayesian Analysis*,
16, 667–718. https://doi.org/10.1214/20-BA1221