Evaluating multi-season occupancy models with autocorrelation fitted to heterogeneous datasets

Technical Summary: Evaluating Multi-Season Occupancy Models with Autocorrelation on Heterogeneous Datasets

Problem Statement
Occupancy models accounting for imperfect detection are standard in ecology for predicting species distributions. However, large-scale naturalist databases (e.g., citizen science) are often highly heterogeneous, characterized by single-visit data, missing observations (spatial gaps), and skewed distributions of survey effort per site. Traditional occupancy models require robust sampling designs (repeated visits within primary occasions) to ensure parameter identifiability (separating occupancy $\psi$ from detection probability $p$). While recent approaches propose using spatial and temporal random effects to represent autocorrelation as a substitute for strict replication ("fractional replication"), it remains unclear how these models perform when applied to real-world datasets dominated by missing data, covariate overlap between occupancy and detection submodels, and spatiotemporal clustering of observations.

Methodology
The authors evaluated a Bayesian multi-season occupancy model with spatial and temporal random effects (based on Doser & Stoudt, 2024) using three simulation studies and empirical data analysis.

1. Simulation Design:

   • Study 1: Replicated the baseline design of Doser & Stoudt (2024) but introduced stronger spatial autocorrelation (slower spatial decay) to test if overestimation of the spatial decay parameter ($\phi$) affects occupancy estimates.
   • Study 2: Shifted from a Bernoulli-distributed survey design to a Poisson-distributed design to mimic the skewed distribution of surveys (including zero surveys) found in butterfly datasets. This study also introduced covariate overlap:
     • Scenario 2-1: Latitude ($L_i$) as an occupancy covariate.
     • Scenario 2-2: Complete overlap of $L_i$ in both occupancy and detection models.
     • Scenario 2-3: Partial overlap of $L_i$.
   • Study 3: Introduced spatiotemporal clustering to simulate phenology (temporal peaks) and observer preferences (spatial "spots" of high effort). This included scenarios with missing data (gaps) both within and outside these clusters.
   • Assessment: The authors assessed parameter identifiability by analyzing the posterior distribution means across 14,400 simulated datasets. Metrics included bias in occupancy estimates ($\hat{\psi}_{it}$), Mean Squared Error (MSE), and contour plots of parameter combinations (e.g., $\hat{\phi}$ vs. $\hat{\sigma}^2$) to detect local identifiability or elongated density shapes. Model misspecification was also tested by simulating data with a probit link and fitting with a logit link.

2. Empirical Data Analysis:

   • The model was fitted to a heterogeneous dataset of 298,389 butterfly records from the Nouvelle-Aquitaine region (France, 2000–2023), aggregated into 1 km² grid cells.
   • Two focal species were analyzed: Polyommatus icarus (Common Blue) and Lycaena dispar (Large Copper).
   • Analyses were conducted on a subset ("Greater Bordeaux area") and the full region to test performance under varying data densities and gap structures.
   • Prior-posterior overlap (PPO) was calculated to evaluate whether data informed the posterior or if priors dominated.

Key Results

• Robustness to Heterogeneity and Overlap: The model demonstrated robustness to skewed distributions of survey numbers (Poisson design) and covariate overlap between occupancy and detection submodels. In scenarios without severe data gaps, the model successfully differentiated the effects of covariates on occupancy and detection, and occupancy estimates remained unbiased relative to the truth.
• Impact of Spatiotemporal Gaps: When spatiotemporal gaps were introduced (Study 3, Scenarios 3-2 and 3-3), the model exhibited significant bias. Site occupancy estimates ($\hat{\psi}_{it}$) were biased toward the average occupancy, which itself was overestimated. This "pull toward the mean" occurred even in data-rich areas within the simulation.
• Identifiability Issues:
  • Spatial Parameters: The spatial decay parameter ($\phi$) and spatial variance ($\sigma^2$) showed weak identifiability. In scenarios with high spatial decay and low variance, the joint posterior density was elongated. The model frequently estimated high $\phi$ (implying short autocorrelation ranges), effectively treating spatial random effects as unstructured noise.
  • Temporal Parameters: Temporal autocorrelation ($\rho$) and variance ($\sigma^2_T$) were consistently underestimated, leading to unstructured temporal random effects.
  • Consequence: These identifiability issues prevented the model from correcting the influence of data gaps. Instead of borrowing strength from neighbors, the model produced nearly constant predictions for unsampled sites, reflecting the biased average occupancy.
• Empirical Findings: In the empirical analysis, the model produced sensible average occupancy trends for the Greater Bordeaux subset but struggled with the full dataset. Parameter estimates for spatial and temporal autocorrelation showed high prior-posterior overlap, indicating the data were insufficient to update the priors. The spatial random effects indicated short-range autocorrelation (or none), and predictions for unsampled areas were flat.

Significance and Claims
The paper claims that while multi-season occupancy models with spatial and temporal autocorrelation are robust to heterogeneous sampling efforts and covariate overlap, they face practical identifiability issues when applied to datasets with severe spatiotemporal gaps.

The authors conclude that:

1. The "fractional replication" approach (using autocorrelation to substitute for repeated visits) works well when data gaps are minimal and the sampling design is not highly clustered.
2. However, in the presence of severe data gaps and spatiotemporal clustering (common in opportunistic datasets), the model fails to correctly estimate variance and autocorrelation parameters.
3. These estimation failures contaminate predictions, causing site-level occupancy to regress toward an overestimated global average, even in areas with data.
4. The authors do not claim the model is useless but caution that its application requires careful assessment of data gaps and parameter identifiability. They note that using informative priors or alternative model structures (e.g., full Gaussian fields via INLA) did not fully resolve the issues in their specific context, suggesting that fundamental limitations exist in estimating these parameters from sparse, heterogeneous data.