Data-augmented models in flocker

When modeling multiple ecologically comparable species simultaneously, occupancy models often assume conditional exchangeability across species and fit species-specific terms as random effects. Data-augmented multi-species models leverage this exchangeability assumption to estimate the number and prevalence of never-detected species within the study region. To do so, the dataset is augmented with a large number of pseudospecies with all-zero detection histories, each species (both the observed species and the augmented pseudospecies) is ascribed a common parameter $\omega$ giving the Bernoulli probability that a given (pseudo)species occurs in the study area, and random-effects exchangeability assumptions are assumed to hold for all species-specific modeled terms (i.e. intercepts and slopes for occupancy and detection).

For a data-augmented multi-species model, we marginalize over the occupancy status of a closure-unit as for a single-season model yielding the unit-wise likelihood $\mathcal{L}_i$ , and we additionally marginalize over the availability of each (pseudo)species, yielding the species-wise likelihood $\mathcal{N}_s = \begin{cases} B(1 | \omega)\prod\limits_{i \textrm{ in } I_s}{\mathcal{L}_i} & \text{if $r_s = 1$} \\ B(0 | \omega) + B(1 | \omega)\prod\limits_{i \textrm{ in } I_s}{\mathcal{L}_i}& \text{if $r_s = 0$} \end{cases}$ where $s$ indexes the species, $\omega$ is the fitted availability probability of a species, $I_s$ is the set of all closure-unit indices $i$ pertaining to species $s$ , $r_s$ is an indicator that takes the value $1$ if there is at least one positive detection of the (pseudo)species in the entire dataset and $0$ otherwise.

Fitting the data-augmented model in flocker requires passing the observed data as a three-dimensional array with sites along the first dimension, visits along the second, and species along the third. Additionally, we must supply the n_aug argument to make_flocker_data(), specifying how many all-zero pseudospecies to augment the data with.

set.seed(2)
library(flocker)
d <- simulate_flocker_data(
    augmented = TRUE
    )
fd = make_flocker_data(
  d$obs, d$unit_covs, d$event_covs,
  type = "augmented", n_aug = 100,
  quiet = TRUE
  )
fm <- flock(
  f_occ = ~ (1 | ff_species),
  f_det = ~ uc1 + ec1 + (1 + uc1 + ec1 | ff_species),
  flocker_data = fd,
  augmented = TRUE,
  cores = 4,
  refresh = 0
  )
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess

Here, the random effect of species is specified using the special grouping keyword ff_species (names beginning with ff_ are reserved in flocker and are not allowed as names for user-supplied covariates).

summary(fm)
#>  Family: occupancy_augmented 
#>   Links: mu = identity; occ = identity; Omega = identity 
#> Formula: ff_y | vint(ff_n_unit, ff_n_rep, ff_Q, ff_n_sp, ff_superQ, ff_species, ff_rep_index1, ff_rep_index2, ff_rep_index3, ff_rep_index4) ~ uc1 + ec1 + (1 + uc1 + ec1 | ff_species) 
#>          occ ~ (1 | ff_species)
#>          Omega ~ 1
#>    Data: data (Number of observations: 25600) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Multilevel Hyperparameters:
#> ~ff_species (Number of levels: 128) 
#>                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept)          1.58      0.27     1.15     2.19 1.00      961     1454
#> sd(uc1)                0.25      0.11     0.06     0.49 1.00     1244     1032
#> sd(ec1)                0.39      0.09     0.23     0.58 1.00     2281     2648
#> sd(occ_Intercept)      1.60      0.33     1.10     2.38 1.00      912     1233
#> cor(Intercept,uc1)     0.74      0.22     0.15     0.98 1.00     2157     2522
#> cor(Intercept,ec1)     0.56      0.20     0.10     0.86 1.00     2935     2678
#> cor(uc1,ec1)           0.57      0.27    -0.07     0.94 1.00      793     1486
#> 
#> Regression Coefficients:
#>                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept           0.48      0.31    -0.14     1.10 1.01      351      528
#> occ_Intercept       0.24      0.34    -0.47     0.89 1.01      525     1094
#> Omega_Intercept    -1.27      0.22    -1.70    -0.85 1.00     7528     3124
#> uc1                 0.01      0.07    -0.14     0.16 1.00      797     1472
#> ec1                 0.07      0.09    -0.10     0.25 1.00     1082     1803
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

flocker enables users to fit data-augmented models using arbitrary brms formulas for the occupancy and detection components. However, we caution that continuous covariates in the occupancy sub-model can lead to pitfalls in interpretation. A seemingly straightforward application, for example, might be to ask how many species are present along an elevational gradient, fitting species-specific quadratic elevation-occupancy relationships. In our experience, a data-augmented model that includes quadratic elevation terms will place an arbitrarily large number of pseudo-species along the gradient, but will do so by placing pseudospecies’ estimated elevational ranges entirely outside the range covered by the sampling effort (Socolar et al 2022). The model is in effect trying to estimate how many species occur in a landscape with elevations ranging from negative to positive infinity. This extrapolation is unprincipled, most obviously because it does not account for hard limits imposed by the physical termina of the gradient (valley floor and mountain peak). Thus, although flocker provides functionality to fit continuous covariates in the occupancy term, we recommend extreme caution in interpreting patterns estimated for never-observed species.

Jacob Socolar

2026-02-24