Prior speci ﬁ cation in Bayesian occupancy modelling improves analysis of species occurrence data

Multi-species biodiversity indicators are increasingly used to assess progress towards the 2020 ‘ Aichi ’ targets of the Convention on Biological Diversity. However, most multi-species indicators are biased towards a few well-studied taxa for which suitable abundance data are available. Consequently, many taxonomic groups are poorly represented in current measures of biodiversity change, particularly invertebrates. Alternative data sources, including opportunistic occurrence data, when analysed appropriately, can provide robust estimates of occur- rence over time and increase the taxonomic coverage of such measures of population change. Occupancy modelling has been shown to produce robust estimates of species occurrence and trends through time. So far, this approach has concentrated on well-recorded taxa and performs poorly where recording intensity is low. Here, we show that the use of weakly informative priors in a Bayesian occupancy model framework greatly improves the precision of occurrence estimates associated with current model formulations when analysing low-intensity occurrence data, although estimated trends can be sensitive to the choice of prior when data are extremely sparse at either end of the recording period. Speci ﬁ cally, three variations of a Bayesian occupancy model, each with a di ﬀ erent focus on information sharing among years, were compared using British ant data from the Bees, Wasps and Ants Recording Society and tested in a simulation experiment. Overall, the random walk model, which allows the sharing of information between the current and previous year, showed improved precision and low bias when estimating species occurrence and trends. The use of the model formulation described here will enable a greater range of datasets to be analysed, covering more taxa, which will signi ﬁ cantly increase taxo- nomic representation of measures of biodiversity change.


Introduction
Targets to stem the loss of biodiversity have been in place globally since 2002 when the Convention on Biological Diversity (CBD) agreed the goal for signatory parties to "significantly reduce the rate of biodiversity loss by 2010". The recognised failure to meet this target was followed by the development of the 'Aichi' targets for 2020 (Convention on Biological Diversity, 2010). The new targets focussed on different facets of biodiversity loss, both direct and indirect, such as awareness of biodiversity, the causes of loss, sustainable land management, the pressures on biodiversity, and the benefits gained from it. To monitor progress towards these goals, a set of biodiversity indicators have been developed to track change in measures related, either directly or indirectly, to these elements (Butchart et al., 2010;Tittensor et al., 2014).
Biodiversity research has, therefore, increasingly focussed on the development of tools to produce robust measures of biodiversity change to accurately measure progress towards these targets (Buckland et al., 2005;Gregory et al., 2005).
Species based indicators are the primary means of monitoring change in the state of biodiversity over time. Several indicators of population change have been developed, at various scales and taxonomic coverage. The Living Planet Index (LPI) is a multi-species indicator that was developed to monitor change in vertebrate abundance at a global scale (Collen et al., 2009). Other examples include the wild bird indicator for Europe (Gregory et al., 2005) and the recent development of butterfly indicators for the UK and Europe (Brereton et al., 2010;Van Swaay et al., 2015). The lack of taxonomic representation is primarily due to dependence on the availability of abundance data from large-scale structured monitoring schemes. In Europe and North America this is limited to birds, mammals, butterflies and moths: elsewhere such schemes are rare. Consequently, current biodiversity indicators are taxonomically biased towards these groups and their ability to act as surrogates of wider biodiversity has been questioned (Rodrigues and Brooks, 2007;Westgate et al., 2014) but rarely evaluated. The lack of taxonomic representativeness is an ongoing problem. If the goal is to understand how biodiversity is changing as a whole, it is important that all groups are represented where possible when such metrics are produced.
One way to achieve greater representation is to use occurrence data for those taxonomic groups that lack abundance data. Occurrence data are "presence-only" records of a species at a known time and location. By also using occurrence data to measure change in biodiversity, it is possible to broaden taxonomic coverage of biodiversity metrics and improve understanding of biodiversity change. For example, the "Dutch LPI" (van Strien et al., 2016) utilises distributional data on dragonflies, fish, mammals, amphibians and butterfly species alongside abundance data. The UK Priority Species Indicator (PSI) uses occurrence data to assess species status for whom abundance data are not available (Eaton et al., 2015;Outhwaite et al., 2015). These indicators have taken advantage of the development of occupancy modelling to incorporate species occurrence data into their assessments (Isaac et al., 2014;Kéry et al., 2010;van Strien et al., 2013).
To date, most applications of occupancy modelling to occurrence data has been limited to well-recorded taxa, such as butterflies (van Strien et al., 2013), dragonflies Termaat et al., 2015;van Strien et al., 2010) and birds (Kamp et al., 2016;Kéry et al., 2010). An exception was the 2015 UK PSI, which used a Bayesian occupancy model framework to analyse a range of taxonomic groups (hymenoptera, bryophytes, carabids, odonata, fish, moths, orthoptera and soldierflies), few of which can be considered to be well-recorded (Isaac and Pocock, 2015). Without occurrence data as an alternative data source, no information would be available for the vast majority of the UK "priority species". The specific occupancy model used for the PSI was that tested by Isaac et al. (2014): they reported high power for estimating species trends compared with alternative methods, and low type I error rates. Isaac et al. explored "high", "medium" and "low" levels of recording intensity, benchmarked against UK and Dutch occurrence datasets. However, occupancy model outputs were useable for only 20% of UK priority species, as most are in taxonomic groups with low levels of recording . Although occurrence data is a vast resource, particularly in Europe and North America, the availability and coverage varies hugely (Meyer et al., 2015). To date, no research has examined the formulation of occupancy models for use with low recording intensity data. To make the most of the occurrence data available the modelling techniques used must be appropriate for the data available and produce outputs with a high precision where possible. The use of current methods on low-intensity data, as shown by the UK PSI, has revealed the need for improvement if they are to be more widely applicable.
One reason for the restricted applicability of current occupancy models is a lack of realistic year-on-year variation in their representation of species occupancy. For example, the model formulation of Isaac et al. (2014) specifies that the occupancy probability of a site is independent from one year to the next. In reality however, for many species the occupancy probabilities in successive years will tend to be similar, with the degree of similarity depending on the species' ecology. This insight can be exploited to constrain the results of an occupancy analysis in a principled manner, with the constraints providing crucial additional "information" that extends the applicability of such techniques to much sparser data sets than has previously been possible. In practice, the constraints are specified using carefully constructed prior distributions within a Bayesian framework. The use of informative and biologically plausible prior distributions can increase confidence in estimates produced from ecological studies (McCarthy and Masters, 2005), although care is required to ensure that the priors do not influence the results unduly. This approach has not been used previously in occupancy trend estimation.
Here, we use the occupancy model framework tested by Isaac et al. (2014) as a base against which to compare alternative specifications that differ in how information about the occupancy state is shared among years. The aim is to determine whether the additional information from sharing across years can advance current modelling practice to improve precision and reduce bias of trend estimates from datasets with low recording intensity. Specifically, we ask: (1) can alternative prior formulations in a Bayesian occupancy model framework improve the precision of species annual occurrence estimates? (2) Do these alternative formulations increase the precision and reduce the bias of species trend estimates compared to the original method tested by Isaac et al. (2014)? If a more appropriate formulation of this occupancy based method can be determined, it will extend the range of taxonomic groups to which occupancy models can be reasonably applied, thus contributing to broadening knowledge on biodiversity status.

Methods
Occupancy models are designed for the analysis of 'presence-absence' data from a collection of sites over time: an occupancy dataset for a particular species typically consists of a set of binary values {y itv } say, where y itv takes the value 1 if the focal species was observed at visit v to site i in year t and 0 otherwise. These elements may be supplemented by other variables, such as sampling effort, or weather, that are potentially related to the probability of observing a species if it is present. To determine whether the use of occupancy models could be improved for the analysis of low recording intensity occurrence data we compared two variants of the Bayesian modelling framework tested by Isaac et al. (2014) to the original model formulation used by those authors (hereafter the 'base model'). We compared model variants using data for ants in Great Britain, for whom the data available is similar to the low recording intensity simulated by Isaac et al. (2014). We also tested the model variants in a simulation experiment to compare their performance with respect to the precision and bias of trend estimates.

The base model
The base model of Isaac et al. (2014) is split into two distinct sub models: a state model and an observation model. The "closure period" (the temporal precision of the state model) is one year; the observation model estimates the probability of detection based on repeat visits within years. The per visit detection history of each species is inferred from records of other species in the assemblage.
The state model, as defined by Eqs. (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. Let ψ it denote the probability that the site is occupied. Then z it has a Bernoulli distribution: Then occupancy probability varies with site and year: where b t and u i are referred to as a 'year effect' and 'site effect' respectively (more details in Eqs. (5) and (6), below). Next, the observation submodel describes how the data were generated. Let p itv denote 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 itself a Bernoulli variable, with conditional distribution modelled as: According to this model, species may only be detected if they are present at a site, so that if z it equals zero then y itv will always be zero. As such, the model assumes there are no false positives (i.e. misidentifications) of the focal species within the data. If a site is occupied (i.e. z it = 1) then (3) gives y itv ∼ Bernoulli(p itv ).
Variation in detection probabilities, p itv , per site, per year and between visits is then modelled as follows: where a t is a year effect and L itv is the list length, defined as the number of species recorded at 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 (Szabo et al., 2010), and it is assumed that the detectability of a species will covary with sampling effort (i.e. generally c > 0). As this model is fitted using Bayesian inference, each of the unknown parameters to be estimated must be assigned a prior distribution describing our knowledge of the system before the data were collected. As is often the case when occupancy models are used to analyse ecological data, the priors of the base model are vague and uninformative, in an attempt to represent a complete lack of information on the system. Isaac et al. (2014) used the following prior formulations: The priors on the variances, σ u 2 and σ a 2 , are specified in the model code via their inverses τ u and τ a . This is standard practice in Bayesian inference, as it is mathematically and computationally convenient. The quantity = τ u σ 1 u 2 is referred to as a precision parameter. The model can be fitted to the data using Markov Chain Monte Carlo (MCMC) techniques, which are designed to produce samples from a posterior distribution of any unknown quantity. Estimates of species occurrence are then derived from the posterior distribution of the parameter z it as the expected proportion of occupied sites for each year.

Problems with the base model
One feature of the model formulation (1)-(10) is that the year effect, b t , is independent between years. This parameter is treated differently from u i and a t , which can be considered as random effects in that they both have a hierarchical prior structure with variance parameters σ u 2 and σ a 2 that are themselves assigned distributions. The variance parameters are estimated along with all other parameters within the model: this in turn enables the site and year effects, u i and a t , to adapt to variation within the data. Since b t does not have a prior distribution including these variance parameters, it does not have this associated flexibility. One implication of this is that it may be necessary to acquire a relatively large number of samples before meaningful estimates can be obtained. Consequently, when there are few data, values drawn from the prior distribution for b t (Eq. (5)) and converted to the probability scale will more than half of the time be either greater than 0.99 or less than 0.01 and, critically, will be independent from year to year. This means that, according to the prior, occupancy state of a site could often switch between very likely to be occupied to highly unlikely to be occupied between years. This does not make biological sense, and the effect is that occupancy estimates can show large fluctuations from year to year when applying these models to datasets of low recording intensity. For example, Fig. 1 shows the base model outputs of annual occurrence for the yellow meadow ant (Lasius flavus), in Great Britain (GB) where estimated occupancy varies by as much as 40% from one year to the next, which is ecologically implausible. Another issue evident in Fig. 1, is the low precision of estimates (credible intervals spanning most of the possible range) in some years. This is because the data are not sufficiently informative to overcome the unrealistic structures implied by the priors. This phenomenon is most pronounced at the start and end of the sequence, since these years tend to have fewer records. Lasius flavus is the third most frequently recorded ant in the dataset of 57 species, with 1862 records, so most models will look worse than this when recording intensity is low. From this model we might tentatively conclude that the species declined during the late 1970s and early 1980s before stabilizing, but the magnitude of change is highly uncertain. For nearly half the years in this sequence, the occupancy estimates are so imprecise as to be uninformative.
Another problem is the presence of boundary effects from the use of restrictive priors on the standard deviation parameters of the hyperpriors (Eqs. (7) and (9)). The standard deviation hyperpriors are set as uniform distributions with limits of 0 and 5. When this parameter is at its maximum value, 5, the associated precision parameter, = τ u σ 1 u 2 , has a minimum value of 0.04. Thus, if the data provide strong information that the true precision is smaller than 0.04, then the posterior distribution cannot adapt fully to this information and instead will be concentrated on the limit set by the prior. Our experience suggests that this problem is not uncommon. Alternative hyperprior options should therefore be considered to address this problem.
We address each of these issues by developing variations on this base model.

Model variant 1: adaptive stationary model
As outlined above, one deficiency with the base model is that the variation in the year effects, b t , is not allowed to adapt to the data. One simple device to allow such adaptation would be to model the prior mean and variance for b t as unknown parameters with hyperprior distributions assigned to them, in the same manner as for the detection sub-model (Eqs. (8) and (9)). In this variant, we replaced Eq. (5) with the following: The priors on other elements of the model remain unchanged. This variant is referred to as the 'adaptive stationary' (AS) model, as it allows the overall mean and variation in year effects to adapt to the data (Gelman et al., 2013, Section 5).

Model variant 2: random walk model
While the adaptive stationary model is expected to be an improvement on the base model, there is still the possible disadvantage that the mean of the year effect is constant over time. Given that the purpose of the analysis is to study long-term changes where these are present, this specification seems potentially limiting for species in rapid increase or decline. As an alternative to the adaptive stationary specification therefore, we consider a variant in which the year effects are allowed to drift systematically over time. To achieve this, we model b t as a random walk (Chandler and Scott, 2011, 5.2, 10.3.1): thus a widely dispersed normal distribution is used for b 1 (representing a priori ignorance about the initial status of a species) and then the subsequent year-on-year changes are modelled using independent normal distributions. This imposes an a priori judgement that the underlying occupancy probabilities are likely to be similar from one year to the next, with the precise degree of similarity controlled by the variance of the year-onyear changes. Random walk priors are commonly used in similar situations (e.g. Chiogna and Gaetan, 2002;Fahrmeir and Lang, 2001;Lee and Shaddick, 2008): they can often improve the precision with which trends can be estimated from sparse data because they allow the sharing of information between time points in a natural way. To implement the approach, we replace equation 5 in the base model with: Normal(0, 100), and (0, 5).
Here the variance σ b 2 controls the extent to which each year's overall mean level is anticipated to deviate from that of the preceding year which ultimately is related to the rate at which a species' distributions increase or decrease. As with the adaptive stationary model, the inclusion of this parameter enables the fitted models to adapt to the data for each species: species that fluctuate little over time will tend to have small values of σ b 2 , while those experiencing substantial variation will have large values. This sharing of information across years will likely increase the precision of occupancy estimates whilst allowing for trends over time. This variant is referred to as the random walk (RW) model. The prior for b 1 has been set with a variance of 10 4 : our experience is that results are generally insensitive to this choice, but that some care is required when data are sparse at either end of the recording period. The relevant issues, and recommendations in the case of sparse data, are discussed in Appendix D.

Hyperprior choice
Both model variants described above have been formulated with the hyperprior choice of the base model, using a uniform distribution with limits of 0 and 5 on the relevant standard deviation parameters. As discussed, this can lead to boundary effects in the precision parameters of the model. As an alternative, we have tested both the adaptive stationary and random walk variants, but replacing these uniform hyperpriors with half-Cauchy hyperpriors. The half-Cauchy is the same as the Student's t distribution with 1 degree of freedom but restricted to positive values, and it has been recommended for use as the default prior for scale parameters within the literature (Gelman, 2006;Polson and Scott, 2012). The distribution is unbounded, unlike the uniform prior, and so does not lead to issues with boundary effects. These differences are tested to determine whether the half-Cauchy is an appropriate replacement for the uniform distribution used on these hyperpriors (Eqs. (12) and (14)) alongside the adaptive and random walk changes. Eqs. (15) and (16) An inverse gamma prior on the variance was also tested, however, estimates frequently failed to converge under this option so it has not been included here.

Model fitting
All models were fitted in a Bayesian mode of inference using JAGS (Plummer, 2009) through R (version 3.1.2, R Core Team, 2017) using the package R2jags (Su and Yajima, 2015). JAGS code for every model can be found in Appendix A.
Each model was fitted for each of the six example species for 50,000 MCMC iterations for three Markov chains with a thinning rate of three. Initial values are set to start the MCMC chains: initial values for z were 1 or 0 depending on whether the species has ever been detected at a site within a year. All other parameter initial values were drawn from a uniform distribution with limits of −2 and 2. These standard settings are not changed within the tests carried out here and are used for each of the alternative model formulations. As the early values of the chain can be highly dependent on these initial values, the first 25,000 iterations are discarded as 'burn-in'.
Convergence of occurrence estimates was determined using the Gelman-Rubin statistic (Rhat) which compares within-chain variance to between-chain variance (Gelman and Rubin, 1992). Convergence was deemed to be acceptable when the Rhat value was below 1.1 (Kéry and Schaub, 2012).
2.6. The data 2.6.1. Species occurrence data Ant data were supplied by the Bees, Wasps and Ants Recording Society (BWARS). This represents one of the many sources of occurrence data available within the UK that are collected opportunistically by volunteer led recording schemes . The ants have relatively few records per species per year (Isaac and Pocock, 2015) and the overall recording intensity (as measured by the distribution of repeat visits across sites) is close to the "low intensity" scenario simulated by Isaac et al. (2014). This made it an ideal taxonomic group for testing how the different model variants would fare for species with few records. Six ant species were selected from the 57 represented within the dataset to demonstrate a variety of data availability per species and varying coverage at a country level. Fig. 2 shows sites where the species was recorded at least once for the period 1970-2013. Records of species within the ant dataset ranged from 2 to 2029 records for the time period assessed, with a median of 107 records over this period.

Ant data preparation
Species records had to meet a number of precision criteria before analysis: only records from 1970 onwards were used, the date of the record had to be known and the location had to be recorded at 1 km 2 precision or better. Data for all species of ant were used: detections of the six focal species were taken directly from the records; non-detections were inferred from the detection of other ant species during visits where the focal species was not recorded. This enabled the generation of detection histories allowing the use of the occupancy models with what would otherwise be presence only data . In order to control for uneven sampling of sites in space and time, sites are filtered based on the number of years with data. Previous studies have used a threshold of three years of data per site Powney et al., 2015). However, Kamp et al. (2016), in an analysis using Danish bird data, reported that precision was lower when filtering on three years compared with unfiltered data. Therefore, we set a threshold of two years thus maximizing the signal in the data whilst excluding sites that contain no information on change.
The final ant dataset consisted of 32,422 records from 18,649 visits to 8,256 unique 1 km grid cells, covering 57 species over 44 years , Table 1. We fitted each of the three model variants, plus those with the alternative hyperprior, for each of the six test species to explore the issues identified above, namely year-to-year fluctuations, the precision of estimates and the presence of boundary effects.

Simulations experiment
We simulated recording data with a known decline to test the three model variants' capacity to estimate known trends. Datasets were generated using the procedure employed by Isaac et al. (2014) using their code with minor changes described below. To ensure that this experiment provided information on how these models might perform in practical situations where the data generating mechanism is unknown, the trends used in generating the datasets do not correspond directly to any of the models being tested.
First, a species occurrence matrix (the 'true' occurrence state at each site) was generated for 25 species (one focal species and 24 non-focal species) across 1000 sites. Species were randomly distributed across sites by sampling 1000 times from a binomial distribution with a species-specific probability of being occupied. The occupancy probability

Table 1
Information on the size of datasets within this study. For the medium and low recording intensity datasets values are means from across the 1500 simulated datasets generated for that recording scenario (500 for each level of change). "Per species" values are calculated across all species in the dataset. "Focal species" refers to the species that were modelled: for the ants these are the six example species. for the focal species in the first year was set to 0.2 (compared with 0.5 in Isaac et al., 2014); probabilities for non-focal species were drawn from a beta distribution with shape parameters of 2 and 2. Within the simulation of these datasets, the focal species was subjected to a set percentage change in probability of occurrence across the 40-year period: we simulated declines rather than increases since applications of such models tend to focus on conservation status for which declines are the primary cause for concern. We set declines of 10, 30 and 50%. The occurrence matrix was then subjected to the "control" recording scenario of Isaac et al. (2014) to generate a dataset representative of species observations. This scenario represents random sampling by a team of virtual observers visiting a certain number of sites each year. The number of visits in a year to each site ranges from 0 to 10, with the probability of receiving n visits given by α·n −2 for 1 ≤ n ≤ 10, so that α is the probability of a site receiving a single visit and the probability of a site not being visited is − ∑ = − α n 1 n 1 10 2 . Levels of α were set to 0.05 and 0.07 to represent low and medium recording intensities (Isaac et al. estimated α = 0.042 for ants and α = 0.07 for dragonflies). After determining the number of sites to be visited in a year, the identities of these sites were selected at random. Visits were then distributed among sites according to the site species richness so that those sites with the most species would receive a greater number of visits. This is representative of real datasets where visits tend to be clustered in biodiverse sites. Each species had a fixed detection probability. That of the focal species was set to 0.2 (compared with 0.5 in Isaac et al. (2014)), those for the non-focal species ranged from 0.16 to 0.88 (for detail on how these values were obtained see Isaac et al., 2014). 500 datasets were simulated for each scenario: the two recording intensities (low or medium) and for each level of decline of the focal species (10, 30 or 50%). These parameters were chosen in order that the total number of records and the number of records of the focal species are both comparable with the ant example data (Table 1).
To each of the simulated datasets (500 for each recording intensity and decline combination), we fitted each of the three model options: the base, the adaptive stationary and the random walk models. The control recording scenario protocol used in our simulations produced datasets in which the records are distributed more evenly among sites, years and visits than in the ants, so models were run for just 20,000 MCMC iterations with the first 10,000 discarded as burn-in. The estimated trend across the 40-year period was calculated from the posterior distribution as the mean difference in estimated occupancy between the first and last year, expressed as a proportion of the first year's occupancy; the precision of this estimate was assessed using 95% credible intervals obtained from the MCMC samples.
To test the performance of the three models in detecting trends in species occurrence, the estimated trend was compared to the proportional change that was used to generate the focal species data (i.e. a proportional decline of either 0.1, 0.3 or 0.5).
Precision of the trend estimates was assessed via the coverage and mean width of the credible intervals of the trend estimates. The coverage was calculated as the proportion of simulated datasets where the true decline falls within the 95% credible intervals of the estimated decline: if the intervals are accurate, the coverage should be around 0.95. The width of the credible interval was calculated for each dataset as the difference between the upper and lower limits of the credible interval, and a mean taken over the 500 datasets. This measure demonstrates how wide the credible intervals are on average, with smaller values indicating narrower intervals i.e. greater precision. An ideal model will yield narrow intervals with a coverage of around 0.95.
Values of these and additional performance measures can be found in Appendix B.

Occupancy estimates of six ant species
The two major issues with estimates from the base model, as illustrated by the Lasius flavus example in Fig. 1, are excessive year-to-year fluctuations and high uncertainty of estimates, particularly at the start and end of the time series. The two alternative model variants are  1970 1980 1990 2000 2010 1970 1980 1990 2000 2010 1970 1980 1990 2000 2010 1970 1980 1990 2000 2010 1970 1980 1990 2000 2010 1970 1980 1990 2000 improved in these respects. Firstly, interannual variation is reduced in both variants, but especially the random walk model (Fig. 3). Secondly, by sharing information across years, the alternative models are able to estimate occurrence with much greater precision. This is particularly evident in years at the start and end of the time period. Estimates from the adaptive stationary model are shrunk towards an overall mean, reducing interannual fluctuations. The random walk model, on the other hand, shares information between neighbouring years resulting in greater precision and smoothed estimates. As a result, trends that were highly uncertain in the base model are more clearly expressed. There are also fewer instances of non-convergence of estimates with the adaptive stationary model and random walk variants than with the base model. Improvements are also evident for the rare and small ranged species, Formica picea and F. aquilonia, where both alternative variants are confident that the species are rare in the early years, unlike the base model where the very wide credible intervals show high levels of uncertainty. Sharing information among years, particularly in the random walk formulation, allows occupancy to be estimated with precision even when there are few records.

Effect of hyperprior choice
The boundary effect problem is clearly evident for the results using the half-uniform hyperprior if we look at the mean estimates of the τ u parameter, Table 2. When this hyperprior is used, the parameter estimate cannot be less than 0.04. Some of the species analysed yield posterior estimates very close to this boundary (see Table 2). This problem disappears when the half-Cauchy hyperprior is used: in most cases when the posterior mean was close to 0.04 under the uniform hyperprior, the half-Cauchy choice yields values that are considerably smaller.
Despite this problem, there were no obvious differences between the mean posterior occupancy estimates produced from the set of models tested when using either the uniform or half-Cauchy hyperpriors, Fig. 4. This suggests that the choice of hyperprior does not have an appreciable influence on the posterior distributions of interest, however some differences in convergence were observed under the adaptive stationary model.  Fig. 4. Plots of species occurrence estimates over time after 50,000 iterations comparing variations of hyperpriors on the standard deviation parameters of the model for six species of ant. Point estimates are coloured according to their Rhat value: estimates with a Rhat below 1.1 (converged) in blue and above 1.1 (not converged) in red. 95% credible intervals are shown in grey. Species detections are plotted along the top of the plots to show data levels over time. Species are ordered by number of records, from highest to lowest. +hC denotes those outputs using the half-Cauchy hyperpriors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Precision and bias of trend estimates from simulated data
In reporting the results of our simulation experiments, we focus on those for which half-Cauchy hyperpriors were used. This is due to their avoidance of boundary effects as discussed earlier. Simulations using half-uniform hyperpriors have also been carried out; the results are presented in Appendix C. Fig. 5 shows the distribution of the 500 trend estimates from each model, for each combination of recording intensity and change scenario. The base model has very low precision, shown by the extremely wide spread of the distribution across the range of possible values. So, individual trend estimates have the potential to be highly under-and over-estimated. The adaptive stationary model produced trend estimates that are the most precise, being concentrated in a small area of the possible outcomes. However, they also show the greatest level of bias, with the distribution of estimates being furthest from the true trend. The random walk model has reduced bias compared to the adaptive stationary and base models and has a greater precision than the base model, (Fig. 5 and Table B1). Although it has a tendency to underestimate change, shown by the distribution being on the more positive side of the true value (dashed line in Fig. 5), it performs better than each of the other model variants.
The adaptive stationary model shows the lowest coverage, with values as low as 0.45 ( Fig. 6 and Table B2). This shows that this model often produces estimates where the 95% credible intervals do not encapsulate the true trend. This is not surprising given that the prior of this model corresponds to an absence of trend. The base model shows good coverage, however the credible intervals are extremely wide (Fig. 6), especially for the low recording intensity. The random walk model performs well on both metrics. It provides precise estimates as shown by the small width of the credible intervals, and coverage is mostly greater than the nominal 0.95 (ranging between 0.916 and 0.998, Table B2).

Discussion
As the deadline for the 2020 biodiversity Aichi targets draws near, it is important that measures of biodiversity change are as representative as possible. Current measures of species status are biased towards certain taxa including birds, mammals and butterflies (Brereton et al., 2010;Collen et al., 2009;Gregory and van Strien, 2010) due to the limited availability of abundance data. Attempts to broaden coverage using occurrence data have enabled the inclusion of other taxa into biodiversity indicators van Strien et al., 2016), however, these studies are still limited to well recorded groups. This study has explored how improvements in modelling methodology can expand the range of taxa for which biodiversity metrics can be obtained, thereby better exploiting the wealth of occurrence data available, increasing the representativeness of the metrics themselves.
Occurrence records offer the opportunity to broaden our understanding of biodiversity change across taxonomic groups if analysed using appropriate statistical methods. Our models of British ants show that it is possible to improve the precision of occurrence estimates gained from low recording intensity data if an appropriate model is used. We have shown how weakly informative priors make it possible to extract useable trend estimates from these methods where the data were previously thought to be limiting. In well-studied countries of Western Europe, this means that multi-species indicators (e.g. Outhwaite et al., 2015;van Strien et al., 2016) can be extended taxonomically to include a broad sweep of invertebrate and plant species. Furthermore, the use of occurrence data could also be used for wellstudied taxa such as dragonflies or butterflies to produce indicators that have a much greater spatial coverage than is currently possible.
Here, we have identified several drawbacks in the base model and proposed alternative model variants that aim to overcome these drawbacks. The first variant, leading to the adaptive stationary model, does not perform well because the model is not able to track genuine trends on the basis of sparse data. This is best visualised by Leptothorax acervorum, where the U-shaped trend is hidden when using the adaptive stationary model; and from the simulation results where trend estimates are strongly biased towards zero. The random walk model, on the other hand, performed well across all criteria with improvements in terms of reduced fluctuations, increased precision and reduced bias of trend estimatesdespite the fact that, as is the case in reality where the true trends are unknown, the structure of the model does not correspond directly to the true trends in the simulation experiment. Of the three model variants considered here therefore, the random walk seems preferable when analysing occurrence data with low recording intensity. In the scenarios considered, the model leads to slightly conservative point estimates of changes in occupancy but provides realistic estimates of uncertainty, in the sense that 95% credible intervals do indeed include the true change in around 95% of cases.
Despite these apparent advantages of the random walk model compared with the others considered here, it should not be used uncritically: in any Bayesian analysis, the implications of prior choices need to be understood. Thus, although our aim was to develop a modelling framework with the flexibility to adapt to a wide range of species, some species dynamics may be harder to capture than others. For example, for species consisting mostly of annual immigrants (e.g. the painted lady Vanessa cardui), occupancy may genuinely fluctuate from year to year: the random walk prior is not specifically designed to represent this behaviour, although it could be approximated by setting a high value for the inter-year variance σ b 2 in Eq. (13) and, as with any Bayesian analysis, the results will become progressively less dependent on the prior with increasing data availability. When analysing individual species or species groups therefore, if data are likely to be very sparse then it is advisable to consider the known dynamics of the species when deciding on the appropriateness of a given prior specification. Nevertheless, when producing large, multi-species analyses such as for biodiversity indicators, a broad approach that works for the majority of species is expedient. A further caveat relates to our claim that the random walk model leads to increased precision of trend estimates. This is evident for the rare ant species F. picea and F. aquilonia (Fig. 3) where the model is highly certain that the species are rare across the time period assessed, including at the start of the time period where data are sparse. This high precision is inferred by the model due to the underlying assumption that year-on-year changes in the log odds for occupancy have a constant variance: this essentially mimics the smoothness of a dynamic occupancy model (van Strien et al., 2013) but with fewer parameters. However, if the inter-year dynamics were to change substantially during periods of sparse data, then the random walk model could not detect this and would therefore produce incorrect assessments of the precision of its occupancy estimates.
The use of random walk priors also has implications for the shape of the distribution for the probability (rather than the log odds) of occupancy, particularly for the initial value in the series. This issue is rather subtle: details, and an exploration of the consequences, are provided in Appendix D. Our investigations suggest that when data are sparse, trend estimates may be sensitive to the precise choice of prior for b 1 (see Eq. (13)); but also that the associated uncertainties tend to be high. To avoid misinterpreting the results of any analysis therefore, uncertainties should always be considered alongside any trend estimates. When data are sparser than for the ant examples considered here, it is possible that both the trend estimates and their uncertainties may be sensitive to prior choices: if this is a concern, the sensitivity should be explored as described in Appendix D.
Another improvement established during this study is the introduction of half-Cauchy hyperpriors to remove boundary effects caused by the use of half-uniform hyperpriors in the base model. The fact that these effects were shown to be present in four of our six example species shows how common a problem this can be. Our investigations suggest that this change in hyperprior formulation makes little difference to occurrence estimates or simulation results, but as a point of principle the half-Cauchy distribution is preferable.
Further improvements could be gained by extending this data sharing approach to the sites as well as the years. This would have potential in improving the precision of estimates even further by sharing information between sites through the use of spatial as well as temporal priors. The downside of this approach would be the additional computational burden: given the vast improvements that have already been demonstrated by changing the year effect prior, this additional burden was not deemed to be worthwhile. Additionally, the inclusion of Julian date into the detection model to describe differences in species detectability throughout the year has been applied elsewhere (van Strien et al., 2013. This was tested in the early stages of model development: however, no positive influence was observed, so the idea was not pursued further.
The improvements detailed in this work focus on the issue of data sparseness when analysing occurrence data. There are further difficulties in assessing species trends using this form of data, which have not been addressed here. Common issues include sensitivity to false positive detections, lack of representativeness of sites surveyed and variation in species detectability between sites. It is possible that low recording intensity data are more susceptible to these problems, which should be considered in future analyses. However, we have been able to show that the alteration of the prior on the year effect of the state model from an independent parameter to a random walk process can greatly improve the estimates of occurrence that can be derived from low-intensity occurrence data. This opens the door for large scale analyses for taxonomic groups that have historically been under-represented in multispecies status studies such as biodiversity indicators.