Reporting a crywolf analysis

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:

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:

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:

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:

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:

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.

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:

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

Here the posterior mean of the occupancy slope moves by at most NA 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 NA, which uses Stan through CmdStan NA. 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(NA, NA) prior truncated at zero; intercepts had Normal(0, NA) priors and slopes Normal(0, NA) priors on the logit scale. We ran NA chains of NA warm-up and NA sampling iterations. The largest R-hat was NA, the smallest bulk effective sample size was NA, and there were NA divergent transitions. A posterior predictive check based on the Freeman-Tukey discrepancy gave a Bayesian p-value of NA.

Occupancy NA with x (slope NA, 95% credible interval NA to NA). The mean false-positive rate per automated occasion was NA (NA to NA). With a narrower gap prior or tighter slope priors, the posterior mean of the slope changed by at most NA.

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