Occupancy models

Occupancy models allow analysis of ‘presence-absence’ data from a set of sites through time. Repeated visits (survey events) within a sampling period (typically a year) are used to estimate the probability that a species is recorded if it is actually present at a site. There are two hierarchically coupled sub-models: one dealing with occupancy (presence vs absence), the other handling probability of observation (i.e. detection or non-detection) on repeat visits within the sampling period. For each species, the probability of detection is inferred from the records of other species in the assemblage. This allows generation of detection histories for use in occupancy models with what would otherwise be presence-only data (Kéry et al., Reference Kéry, Royle, Schmid, Schaub, Volet, Häfliger and Zbinden2010). Detection histories are sequences of presences or absences of a given species for a series of visits during a sampling period.

The occupancy model used here is that available within the Sparta software package (August et al., Reference August, Powney, Outhwaite, Harrower, Hill, Hatfield, Mancini and Isaac2021) developed and provided by the Centre for Ecology and Hydrology (CEH). The original approach developed by Isaac et al. (Reference Isaac, van Strien, August, de Zeeuw and Roy2014) was refined and improved by Outhwaite et al. (Reference Outhwaite, Chandler, Powney, Collen, Gregory and Isaac2018). Description and tests of these improvements and the specific details of the model are provided in Outhwaite et al. (Reference Outhwaite, Chandler, Powney, Collen, Gregory and Isaac2018), but briefly, the improvements included: (i) introduction of a random walk component for the occupancy priors to allow more precise estimates of occurrence and smoother trends, particularly where data are sparse, because information is shared between adjacent time periods in a biologically meaningful way; and (ii) replacement of prior estimates taken from a uniform distribution to those from a half-Cauchy distribution (Outhwaite et al., Reference Outhwaite, Chandler, Powney, Collen, Gregory and Isaac2018).

The base model used in Sparta (Isaac et al., Reference Isaac, van Strien, August, de Zeeuw and Roy2014) is split into two distinct sub-models: a state model and an observation model. The temporal precision of the model specified by the ‘closure period’ is one year. The observation sub-model estimates the probability of detection based on repeat visits to a site within years. The per-visit detection history of each species is inferred from records of other species in the assemblage.

The state model (equations (1) and (2)), describes the true occupancy state, z_it, of site i in year t. This will be 1 if occupied or 0 if unoccupied. ψ _it is the probability that a site is occupied and z_it has a Bernoulli distribution:

(1)$$z_{it}\sim {\rm Bernoulli( } {{\psi}} _{it}) , \;$$

where the logit of the probability of occurrence ψ _it, varies with year and site:

(2)$${\rm logit}( {{\psi}_{\it it}} ) = {\rm log}\left({\displaystyle{{{\psi}_{\it it}} \over {1- {\psi}_{\it it}}}} \right) = b_t + u_i$$

with b_t and u_i denoting year- and site-effects, respectively.

According to the spatial scale of interest, the state model year-effect was partitioned in different ways to allow the estimation of occupancy for different areas within the UK. For instance, the year effect was split into components representing England, Wales and Scotland when considering the scale of ‘country’. Similarly, when considering smaller scales, the year effect was split into areas representing Isles of Scilly, Cornwall, North Devon, South Devon and Dorset. A third division separated visits by designation status, i.e. whether they occurred inside or outside MPAs in SW England (i.e. in the waters of Cornwall, Devon and Dorset). Visits were coded by the different categories by using ‘join by location’ in QGIS with polygons representing the different territorial waters, areal waters or MPA boundaries.

This means that, instead of having a single year effect in the state model as shown in equation (2), there is a year effect associated with each country, area or designation status. In the case of countries, let r(i) be the country (England, Scotland or Wales) in which site i is located, then:

(3)$${\rm logit}\,( {{\psi}} _{it}) = {\rm log}\left({\displaystyle{{{{\psi}}_{it}} \over {1-{{\psi}}_{it}}}} \right) = b_{tr( i) } + u_i$$

where b _tr(i) is the year effect for year t in country r in which site i occurs.

The observation sub-model describes the data collection process and is conditional on the true occupancy state z_it. p_itv denotes the probability that a species will be observed on a single visit, given that it is present at the site (z_it = 1). Then the observation parameter y_itv is described as being drawn from a Bernoulli distribution conditional on the true occupancy state:

(4)$$y_{itv}\vert z_{it\;}\sim {\rm Bernoulli( }p_{itv}.z_{it}) $$

According to this model, species may only be detected if they are present at a site, so if z_it equals zero, then y_itv will always be 0. If a site is occupied (i.e. z_it = 1) then equation (4) gives y _itv ~ Bernoulli(p _itv).

Variation in detection probabilities (p_itv) per site, per year and between visits is then modelled as follows:

(5)$${\rm logit}( {\,p_{itv}} ) = \log \left({\displaystyle{{\,p_{itv}} \over {1-p_{itv}}}} \right) = a_t + c{\rm log}L_{itv}, \;$$

where a_t is a year effect and L_itv is the list length, defined as the number of taxa recorded during a single visit v to site i in year t. Parameter c represents the relationship between overall sampling effort and the detection probability of the focal species. List length is used as a proxy for sampling effort per visit (Franklin, Reference Franklin1999; Szabo et al., Reference Szabo, Vesk, Baxter and Possingham2010), and it is assumed that the detectability of a species will co-vary with sampling effort (i.e. generally c > 0). In a previous application of this approach (Outhwaite et al., Reference Outhwaite, Powney, August, Chandler, Rorke, Pescott, Harvey, Roy, Fox, Roy, Alexander, Ball, Bantock, Barber, Beckmann, Cook, Flanagan, Fowles, Hammon, Harvey, Hepper, Hubble, Kramer, Lee, MacAdam, Morris, Norris, Palmer, Plant, Simkin, Stubbs, Sutton, Telfer, Wallace and Isaac2019), list length was included using a categorical specification, with lists of length 1, 2–3 or 4+. In the present study, a continuous specification was retained due to the typically large list-lengths (often up to 50) relative to samples from the terrestrial environment, where the three categories of short list length seemed inappropriate.

For Bayesian models, each of the unknown parameters to be estimated needs to be assigned a prior distribution that reflects our understanding of the system, in advance of data collection. When analysing ecological data, we often know very little about the system and so vague priors are typically used. With the exception of the year effect (b _t), all parameters were assigned vague priors. The parameters describing the distributions of these priors (i.e. hyperpriors, e.g. σ ²) were assigned half-Cauchy distributions based on Cauchy distributions with location 0 and scale parameter 1 (which is the same as the modulus of Student's t distribution with 1 degree of freedom (d.f.)), as recommended by Outhwaite et al. (Reference Outhwaite, Chandler, Powney, Collen, Gregory and Isaac2018).

(6)$$u_i\sim {\rm Normal}\,( 0, \;\sigma _u^2 ) , \;\quad {\rm where}\;\sigma _u\sim \;\vert {{\rm Student}\ t\;{\rm on}\;1\;{\rm d}.{\rm f}.} \vert $$

(7)$$a_t\sim {\rm Normal}\,( \mu _a, \;\sigma _a^2 ) , \;\quad {\rm where}\;\mu _a\sim {\rm Uniform}\,( 0, \;\,100) $$

and

(8)$$\eqalign{\sigma _a\sim & \vert {{\rm Student}\ t\;{\rm on}\;1\;{\rm d}.{\rm f}.} \vert \cr c\sim & {\rm Uniform}\,(-10, \;\,10) } $$

The year effect, b _t, uses a ‘random walk’ prior, which improves the precision of occupancy estimates by allowing for yearly variation among probability distributions. Thus:

(9)$$b_t\sim \left\{{\matrix{ {{\rm Normal}\,( \mu_{br}, \;\;{10}^4) \quad {\rm for}\;t = 1} \cr {{\rm Normal}\,( b_{t-1r}, \;\sigma_{br}^2 ) \quad {\rm for}\;t > 1} \cr } } \right.$$

where,

(10)$$\mu _{br}\sim {\rm Normal}\,( 0, \;\,100) \;{\rm and}\;\sigma _{br}\sim \vert {{\rm Student}\ t\;{\rm on}\;1\;{\rm d}.{\rm f}.} \vert $$

The Sparta occupancy model selected (using ‘random walk’ design, with half-Cauchy distribution of priors) gives greater precision than other models applied previously (Outhwaite et al., Reference Outhwaite, Chandler, Powney, Collen, Gregory and Isaac2018, Reference Outhwaite, Powney, August, Chandler, Rorke, Pescott, Harvey, Roy, Fox, Roy, Alexander, Ball, Bantock, Barber, Beckmann, Cook, Flanagan, Fowles, Hammon, Harvey, Hepper, Hubble, Kramer, Lee, MacAdam, Morris, Norris, Palmer, Plant, Simkin, Stubbs, Sutton, Telfer, Wallace and Isaac2019). This approach appears to be the most appropriate for estimating population trends, particularly where the number of records is not large (as happens where sampling intensity is not great and/or areas of interest are small) or inhabited range is small. This model also provides estimates of uncertainty (Outhwaite et al., Reference Outhwaite, Chandler, Powney, Collen, Gregory and Isaac2018, Reference Outhwaite, Powney, August, Chandler, Rorke, Pescott, Harvey, Roy, Fox, Roy, Alexander, Ball, Bantock, Barber, Beckmann, Cook, Flanagan, Fowles, Hammon, Harvey, Hepper, Hubble, Kramer, Lee, MacAdam, Morris, Norris, Palmer, Plant, Simkin, Stubbs, Sutton, Telfer, Wallace and Isaac2019). The models were implemented in the R programming environment (version 4.0.3, R Core Team, 2021) and fitted with the JAGS programme for analysis of Bayesian hierarchical models (Plummer, Reference Plummer2003).

One assumption of an occupancy model is that recording one taxon tells you something about a taxon that was not recorded. In terrestrial recording, records are typically made within restricted taxonomic contexts (e.g. dragonflies or water-beetles or lichens). This is not so within Seasearch. During Seasearch surveys, as many organisms as possible that can be reliably identified in situ are recorded, irrespective of taxonomic classification. This provides a ‘complete list’ on each visit (Isaac & Pocock, Reference Isaac and Pocock2015). A ‘complete list’ is where a record is made for each species observed (rather than just incidental records of ‘interesting’ species). This is not the same as saying that all species present and observable were actually observed, but rather just that all species observed were reported. This means that a very large pool of records for many taxa are available to be used by the model when estimating probabilities of detection.

The focus of this study was crawfish, brown crab and lobster. The model was fitted and population trends generated for each of these species with 30,000 Markov chain Monte Carlo (MCMC) iterations for three Markov chains, with a thinning rate of three. The first 10,000 iterations were discarded as ‘burn-in’ since early values of the chain can be influenced strongly by the initial values chosen for parameters (Outhwaite et al., Reference Outhwaite, Chandler, Powney, Collen, Gregory and Isaac2018). Whether or not estimates of occurrence had converged by the end of the iterations was determined using the Gelman-Rubin statistic (Rhat), which compares within-chain variance to between-chain variance (Gelman & Rubin, Reference Gelman and Rubin1992). Convergence was considered acceptable when the Rhat value was below 1.1 (Kéry & Schaub, Reference Kéry and Schaub2012).

This analysis allowed, to my knowledge, the first estimation of species occupancy for any benthic marine species in the UK.

Data collection and processing

The raw data used here were occurrence records collected by Seasearch volunteers, held by Seasearch and available through the NBN Atlas (https://nbnatlas.org/; https://doi.org/10.15468/kywx6m; https://doi.org/10.15468/0hyjxi; https://doi.org/10.15468/4us2hk; https://doi.org/10.15468/pyugge; https://doi.org/10.15468/0ppp4p; https://doi.org/10.15468/mxkbcg). For some taxa, publicly available data are provided at a resolution coarser than the original records used here. Data collection protocols have remained unchanged since the start of this century. This means all taxa can be treated in the same analysis and that Seasearch records from 2000 onwards can be used in their entirety, with consequent greater ability to estimate detectability of the modelled species (than if smaller taxonomic subsets were used).

The dataset was first filtered for desired spatial and temporal coverage and standardized to meet the requirements of the model such that:

• • The location of the record was within the waters of the British Isles and adjacent seas (i.e. including England, Scotland, Wales, Ireland, Isle of Man and the Channel Isles).

• • Only records more recent than 01/01/2000 were used.

• • Taxonomic determinations above the level of Family were eliminated.

• • All records were marked as being ‘alive’ and ‘certain’.

• • The date of the record was known to the day (true for all Seasearch records).

• • The position of the record was known to at least 1 × 1 km precision (true for all Seasearch records).

• • Veracity of taxonomic determinations was reviewed by regional experts as per the Seasearch QC process (Bolton, Reference Bolton2018b).

• • When, within any monospecific genus (according to Marine Species for the British Isles and Adjacent Seas – MSBIAS; www.marinespecies.org/msbias/), there existed some records determined only to genus, all entries were altered to the species.

• • Seasearch surveys often include multiple samples (from different habitats) at a single location. Occupancy models tend not to resolve this level of detail yet, so all records from a site are treated together, with duplicate records from a dive being eliminated (i.e. where the same taxon was recorded in multiple samples). If data are finely enough resolved then it may be possible to fit models that treat sites and habitat independently.

• • Multiple visits to the same site on the same day (i.e. multiple divers submitted separate records from the same 1 × 1 km grid square on the same day) were treated as ‘replicates’ in the model.

Data-fields used were the simple ‘what, when, where’ required by the Sparta model, with ‘what’ being the taxon recorded, ‘when’ the calendar date it was recorded and ‘where’ was the position. Recorded positions are considered accurate to about 70 m, but were scaled-up to 1 × 1 km grid cell resolution in order to maximize the number of visits to a ‘site’ within a sampling period (i.e. year) whilst retaining spatial resolution fine enough to be useful. These 1 × 1 km cells coincided with those of the British National Grid, but were extended offshore where necessary.

Tools within Sparta were then used to create detection histories from the standardized data for each site. Annual estimates of occupancy were extracted separately for countries (England, Wales and Scotland), for areas in the south-west of England (Isles of Scilly, Cornwall, North Devon, South Devon and Dorset) and for areas in south-west England that were or were not designated as MPAs.

Changes in the probabilities of annual occupancy through time define trends in species occupancy, where occupancy is the proportion of surveyed grid-cells where the modelled species is present. Species trends were estimated as the annual percentage growth-rate using the first and last years for which that species had records. As ψ _it is the annual occupancy estimate for a site i in year t, the mean occupancy for all (n) sites in year t is ψ¯ _t = $( \mathop \sum \nolimits_{i = 1}^n$ψ _it/n), the yearly growth rate for a site is λ_it = (ψ _it/ψ _it−1), the average growth-rate for a given year is λ¯_t = (ψ¯ _t/ψ¯ _t−1), the total growth-rate is λ_T = (ψ _f/ψ _s) and the annual growth-rate as a percentage is λ = ((ψ _f/ψ _s)^1/y − 1) × 100, where f was the occupancy in the final year, s was the occupancy in the starting year and y was the number of years. The calculation for annual growth-rate assumes a monotonic change in growth, where occupancy changes by the same proportion each year.

For each country or area, annual growth-rate in species occupancy was calculated for the time-period between the first and last records for each of 1000 samples from the posterior distribution of estimates for annual occupancy. These were then used to calculate mean annual growth-rates, an equal-tailed 95% credible interval (C.I.) and a measure of precision (estimated as the inverse of the variance of the 1000 values for growth-rate). Credible intervals (rather than confidence intervals) are used, as appropriate for Bayesian analyses. Using first and last records for a species rather than always using the full 21 years (i.e. 2000–2020) avoided extrapolation beyond the actual data. The mean estimated occupancies (with 95% C.I.) for each year were then plotted to show the population trend for that species in that area. Samples of the posterior distribution of occupancy for each area per species per year are provided (Supplementary materials 1; one file per species per set of areas).

To help reduce temporal bias in geographic coverage caused by uneven sampling, sites can be filtered according to the number of years when records were made. Such filtering will increase spatial bias, which may or may not be preferable depending on the goal of the modelling. Early studies using this method required sites to have records in at least 3 years (Isaac & Pocock, Reference Isaac and Pocock2015; Powney & Isaac, Reference Powney and Isaac2015). Kamp et al. (Reference Kamp, Oppel, Heldbjerg, Nyegaard and Donald2016) analysing trends in populations of Danish birds demonstrated, however, that precision is greater when using unfiltered data. So, to maximize the extractable information whilst excluding sites with no information about change, a minimal threshold of 2 years of data was set when selecting sites (following Outhwaite et al., Reference Outhwaite, Powney, August, Chandler, Rorke, Pescott, Harvey, Roy, Fox, Roy, Alexander, Ball, Bantock, Barber, Beckmann, Cook, Flanagan, Fowles, Hammon, Harvey, Hepper, Hubble, Kramer, Lee, MacAdam, Morris, Norris, Palmer, Plant, Simkin, Stubbs, Sutton, Telfer, Wallace and Isaac2019).

To help ensure that ‘adequate’ data are used by the model, criteria were set, which if not met, precluded the model from running. In previous studies (Outhwaite et al., Reference Outhwaite, Chandler, Powney, Collen, Gregory and Isaac2018, Reference Outhwaite, Powney, August, Chandler, Rorke, Pescott, Harvey, Roy, Fox, Roy, Alexander, Ball, Bantock, Barber, Beckmann, Cook, Flanagan, Fowles, Hammon, Harvey, Hepper, Hubble, Kramer, Lee, MacAdam, Morris, Norris, Palmer, Plant, Simkin, Stubbs, Sutton, Telfer, Wallace and Isaac2019) these criteria were a maximal gap in records of 10 years and a minimum of 50 records in the area of interest. These may be perceived as rather low thresholds. Although low thresholds provide access to more data, more stringent criteria are used here, i.e. a maximal gap in records of 5 years and a minimum of 100 records in the area of interest. This is justified because smaller amounts of good quality data are better than lots of poor data (Meng, Reference Meng2018) and can offset somewhat the biases prevalent in citizen science data which greatly reduce the effective sample size.

Graphical checks were used to increase confidence in the suitability of the model selected and of the trends produced. These included inspection of:

1. (i) patterns of potential bias in when or where records were collected;

2. (ii) patterns in simulations of probabilities of occupancy;

3. (iii) traceplots of individual chains to check that parameters have converged.

Patterns in data collection can mask, enhance or otherwise confound any observed trends in occupancy. To assess whether different forms of spatial or temporal bias in collection of Seasearch records were explored using the R package occAssess (Boyd et al., Reference Boyd, Powney, Carvell and Pescott2021). Forms of bias considered were variation in: sampling intensity through time, taxonomic coverage, proportion of records made to species, taxonomic bias, spatial coverage and spatial clustering. Simulations of occupancy and traceplots were inspected for a subset of model outputs, specifically those showing notable and interesting trends (details provided in Supplementary materials 2 & 3).

Using a desktop PC with Intel^® Core™ i7-8700 CPU @ 3.20 GHz, 32 GB RAM, 64-bit operating system, ×64-based processor, running Windows 10 Pro, the code to run models for three species partitioned into three different sets of spatial scales (countries, areas or MPAs) took just under 14 h to complete. Supplementary materials 1 contains the raw data, supporting files and code to run the model, plus key outputs from the model (e.g. posterior estimates of occupancy, population growth-rates and Rdata files used for Supplementary materials 3).