Accounting for spatial autocorrelation and environment are important to derive robust bat population trends from citizen science data

Monitoring wildlife populations is essential if global targets to reverse biodiversity declines are to be met. Recent analysis of data from the UK ’ s long-term National Bat Monitoring Programme (NBMP) suggests stable or increasing population trends for many bat species, and these statistics help inform progress towards national biodiversity targets. However, although based on robust citizen science survey designs, it is unknown how sensitive these trends are to spatial and environmental biases. Here we use Bayesian hierarchical modelling with integrated nested Laplace approximation (INLA), to examine the impact of these types of biases on the population trends using relative occupancy of four species monitored by the NBMP Field Survey in Great Britain (GB): Pipistrellus pipistrellus, P. pygmaeus, Nyctalus noctula and Eptesicus serotinus . Where possible, we also disaggregated trends to national levels using the best model per species to determine if national differences in trends remain once sampling biases are accounted for. Although we found evidence of spatial clustering in the NBMP Field Survey locations, the previously reported GB-wide population trends are broadly robust to spatial autocorrelation. In most species, accounting for spatial autocorrelation and species-environment relationships improved model fit. The nationally disaggregated models highlighted that GB-wide trends mask differences between England and Scotland, consistent with previous analysis of these data, as well as illustrating large gaps in survey effort, especially in Wales. We suggest that although bat population trends were found to be broadly robust to sampling biases present in these data, small differences could propagate over time and this impact is likely to be more severe in less structured citizen science data. Therefore, ensuring trends are robust to sampling biases present in citizen science datasets is critical to effective monitoring of progress towards biodiversity targets, managing populations sustainably, and ultimately a reversal of global declines.


Introduction
Despite repeated international commitments to improve the status of biodiversity and ecosystems, average extinction risk continues to rise in many taxa, and the bulk of the Convention of Biological Diversity 2020 Aichi targets have been missed (Convention on Biological Diversity, 2019; Díaz et al., 2019;Mace et al., 2018).More ambitious action is needed, not only to halt biodiversity loss but to restore biodiversity and reverse declines (Mace et al., 2018).Measurable biodiversity indicators have a critical role to play to enable specific goals to be reached within a limited timeframe.Current important global indicator metrics, such as the Living Planet Index (Collen et al., 2009;McRae et al., 2017), depend on the collation of regional, national, and local scale biodiversity monitoring.Birds have been particularly well monitored with structured surveys, such as the North American Breeding Bird Survey (Sauer et al., 2017) and the British Breeding Bird Survey (Freeman et al., 2007).These monitoring data have been important for evidencing the impact of anthropogenic environmental changes, such as agricultural intensification and climate change, on biodiversity (Burns et al., 2016;Eglington & Pearce-Higgins, 2012;McDermott Long et al., 2017).Non-structured datasets compiled of mainly opportunistic occurrence records have also been successfully used to produce trends in occupancy for insects, adding important evidence to the insect decline debate (Outhwaite et al., 2020).
Large-scale, long-term monitoring is often only possible if carried out by volunteer citizen scientists.However, relying on volunteers can generate spatial bias towards areas of higher human population density, more easily accessible sites, habitats favoured by volunteers or those perceived to be more likely to contain the target taxa (Isaac & Pocock, 2015).Stratification of survey location, for instance based on habitat type, partly mediates the issue (Buckland & Johnston, 2017), still, volunteers cannot be compelled to monitor specific sites.An assessment of the data underpinning UK biodiversity indicators reported that survey sites for most schemes are not randomly distributed in space, nor between individual countries in Great Britain (GB) or landcover types (Isaac et al., 2016).As data collected by these surveys are used to inform policy, post-survey statistical techniques that robustly control for spatial biases when modelling population trends are therefore essential (Isaac et al., 2014;Isaac & Pocock, 2015).Generalised additive models (GAMs) have been widely used to model non-linear population trends over time (Barlow et al., 2015;Collen et al., 2009;Eglington & Pearce-Higgins, 2012) due to their implementation of splines that account for temporal autocorrelation, where time points close to each other are likely to be more similar than those further away.However, these methods do not easily control for spatial clustering of survey locations and substantial spatial autocorrelation may remain between data points, meaning data from proximal sites are more likely to be similar.Furthermore, for mobile animals the typical movement range can be larger than the distance between survey sites, so it is possible that individuals from one population might be recorded in multiple study sites leading to an inflation of the true sample size.
Approaches to deal with spatial autocorrelation include thinning the data, geographically weighting certain regions (Kupfer & Farris, 2007) or grouping points close together into regions and fitting these as random effects with separate intercepts.Other approaches include incorporating an effect for preferential sampling where the stratified sampling design was not followed (Conn et al., 2017).More recent approaches use Bayesian techniques such as using Markov chain Monte Carlo (MCMC) or Integrated Nested Laplace Approximation (INLA) to extend generalised additive techniques to include spatial dependency (Fahrmeir et al., 2004;Fahrmeir & Lang, 2001;Rodhouse et al., 2011).Whilst flexible, MCMC methods commonly have a long computational time (Fahrmeir et al., 2004) limiting the feasibility of the approach for large datasets.INLA is a deterministic algorithm for Bayesian inference designed for latent Gaussian models able to accommodate linear and nonlinear effects, including nonlinear temporal trends and spatial random effects.The use of approximation and integration drastically reduces typical Bayesian computational times (Beguin et al., 2012;Rue et al., 2009).To date, the INLA modelling method has not been widely applied to population models.These spatially explicit methods open-up the possibility of accounting more robustly for the spatial structure of sampling sites, resulting in more accurate estimates of population trends.
Bats are considered valuable bioindicators (Jones et al., 2009;Russo & Jones, 2015), and in the UK are monitored by the Bat Conservation Trust's (BCT) National Bat Monitoring Programme (NBMP) in accordance with both national and, until recently, European Union law.The NBMP consists of three volunteer-led survey methods, winter hibernacula counts, summer roost emergence counts, and summer acoustic surveys.Analysis of these data suggest stable or recovering population trends for 11 species, however, rarer and more habitat-specialised bat species are not currently monitored (Bat Conservation Trust, 2021).The NBMP Field Survey is a structured summer acoustic survey and a key contributor to British bat population trend estimates for four species: Pipistrellus pipistrellus, Pipistrellus pygmaeus, Nyctalus noctula and Eptesicus serotinus.Volunteers survey 1 km 2 sites randomly stratified based on the Institute of Terrestrial Ecology landcover classes (Bat Conservation Trust, 2021).GAMs are used to model population trends across the UK and within each of the four countries where data allow (Bat Conservation Trust, 2021).GB data are weighted to allow for the different sampling rates in England, Scotland, and Wales by weighting each site in proportion to the ratio of non-upland area to number of sites surveyed for the relevant country.Weighting is not applied to E. serotinus, as this has a restricted range within Britain.However, the NBMP Field Survey is still affected by spatial biases within countries (Isaac et al., 2016).For example, clusters of surveys are conducted in areas of higher human population densities, such as in the southeast of England, and the current analytical approaches may not adequately remove spatial autocorrelation or account for sampling biases within each nation (Kupfer & Farris, 2007).
Other sources of bias in the data may also be impacting the estimated bat population trends.Factors affecting species-specific detection, such as temperature during the survey, are already included as covariates in the GAMs and a site effect partially controls for inter-site variation in habitat, assuming these remain over time (Barlow et al., 2015).Site occupancy is affected by environmental factors, such as the percentage of woodland and urban landcover, and climatic conditions (Boughey et al., 2011;Dietz et al., 2009;Lintott et al., 2016;Santini et al., 2018).Due to uneven sampling across GB, some landcover classes and climatic conditions are under or over-represented in the NBMP Field Survey data, such as the proportion of urban landcover.If populations are increasing more at sites with a greater proportion of urban landcover the estimated trends will be misrepresentative of the true GB trends.Accounting for the landcover type or configuration, and mean climatic conditions annually is therefore important in ensuring population trend estimates are nationally representative.Spatial autocorrelation is also unlikely to be uniform across GB, with possible regional or local differences in the impact of environmental factors on occupancy.Although population trends are currently reported at UK, GB, and country-level, where data availability allows (Bat Conservation Trust, 2021), assessing whether sampling location and environmental heterogeneity are impacting the trends at national scales has important implications for national conservation policy decision-making.
Here, we use data collected on common bat species P. pipistrellus, P. pygmaeus, N. noctula and E. serotinus by the NBMP Field Survey between 1998 and 2016 in GB to examine whether accounting for spatial autocorrelation, or environmental factors impact the direction or magnitude of the British or nationally disaggregated population trends.We use a non-spatial Bayesian hierarchical framework with INLA to first create a base model for each species using a similar structure to the GAMs reported in Barlow et al. (2015) and then compare these to a series of subsequent models that use the Bayesian framework to account for spatial and environmental biases in the NBMP Field Survey data.Using the best performing model per species, trends are then disaggregated to country-level to determine if national differences in trends remain once biases are accounted for.

Methods
2.1 Data 2.1.1NBMP Field survey data We used data collected by NBMP Field Survey volunteers between 1998 and 2016 (number of surveys = 5163) in GB (Fig. 1).Field Surveys were carried out during July -September using heterodyne bat detectors, which were tuned to specific sound frequencies.Trained volunteers were randomly allocated a 1 km 2 grid square close to where they lived, based on a random stratification of UK landcover classes (Bat Conservation Trust, 2021).Volunteers walked a 3 km transect on two separate nights, separated by at least five days, starting 20 min after sunset.Surveys were not carried out during weather conditions known to adversely affect bat activity, such as temperatures below 7 • C at sunset, heavy rain, and strong winds.The transect was split into 12 parts; volunteers made 12-point counts at intervals along the transect counting the number of P. pipistrellus and P. pygmaeus passes whilst the detector was tuned to 50 kHz as the echolocation call peak frequencies of these species are 45 kHz and 55 kHz, respectively.N. noctula and E. serotinus were monitored during the 12-walks in between point counts whilst the detector was tuned to 25 kHz as the echolocation call peak frequencies of these species are 25 kHz and 25-30 kHz, respectively E. serotinus were not surveyed in Scotland as this is a southerly restricted species in GB (Bat Conservation Trust, 2021).The data were analysed as a series of 12 presence-absences per transect, per species and we assumed no false-negatives or mis-identifications were present.In most cases survey metadata were also recorded by the volunteer, including detector type (Table SI.1), temperature, wind speed and survey duration (Table SI.2).

Environmental data
We collected the following environmental variables for the UK (Table 1) at a yearly resolution to reflect change over time, selected based on prior knowledge of bat species' environmental requirements (Boughey et al., 2011;Catto et al., 1996;Davidson-Watts et al., 2006;Dietz et al., 2009;Lintott et al., 2016): (1) Annual percentage landcover type at 1 km 2 for broadleaved and mixed woodland, needle-leaved woodland, freshwater, grassland, agriculture/cropland and urban were calculated from the 300 m 2 European Space Agency Climate Change Initiative Land Cover (ESA CCI-LC) project (ESA, 2017).The groupings of ESA CCI-LC classes used to aggregate to broad categories used in these analyses is available in Appendix B Table SI.3;(2) Annual mean spring precipitation (March-May) (ml/day); (3) Annual mean summer surface air temperature (March-May) ( • C); (4) Annual mean winter surface air temperature (December -February) ( • C).We used climatic variables from Climate, Hydrology and Ecology research Support System and we calculated annual mean seasonal values from monthly data at a 1 km 2 resolution (Robinson et al., 2017).We calculated the total percentage coverage or overall mean of each environmental variable across GB and at NBMP Field Survey sites highlighting the over-representation of some variables, such as urban landcover, in the data (Table 1).

Statistical analysis
All analyses were carried out in R (R Core Team, 2020).

Impact of spatial and environmental bias on relative occupancy
We conducted Bayesian hierarchical modelling using the R-INLA package (Rue, Martino, & Chopin, 2009;Lindgren & Rue, 2015).Our formulation was broadly similar to a traditional occupancy detection model, (MacKenzie & Bailey, 2004).However, we did not hierarchically model the detection probabilities per site as we were interested in relative change over time, rather than predicting occupancy at nonsurveyed locations.Three models were constructed to estimate an index of relative occupancy using the number of points/walks per survey where a species was recorded as present (presence-absence at points/ walks, n = 12), y pa .A binomial distribution was assumed, and a logit link function applied.The probability of a species' occupancy and the number of points/walks with an occurrence per survey, i, are denoted by p i, and n i , respectively: First, we constructed a model for each species with the same structure as currently used for GAMs to report the population trends, including weighting the observations by the proportion of non-upland area per number of sites per country (Barlow et al., 2015;Bat Conservation Trust, 2021); referred to as the "Base model".Under the Base model p i is modelled as a function of survey specific covariates and random effects, where the model intercept is α 0 ; β is the vector of regression parameters and X i is a matrix of the survey specific covariates, with a vector of linear coefficients (Table SI.2).Continuous variables were standardised to the variation around the mean (mean = 0, sd = 1).A second-order random walk was used as a smooth term for the random effect of count year (t i ).
Site was modelled as an independent and identically distributed random effect (δ i ).We confirmed that the resulting trends from the Base model are broadly consistent with reported GAM trends for each species (Fig. SI.1) (Bat Conservation Trust, 2021).This Base model was used to assess the impact of accounting for spatial bias and then environmental bias on the resulting population trends, overall model fit and predictive accuracy.See Text SI.1 for details of the priors used on random effects.We confirmed spatial autocorrelation occurred between sites using the Moran's Index between sites (p ≤ 0.01, observed = 0.01, expected = − 0.00).We constructed an INLA model (referred to as "Base model + SPDE") that included a spatial random effect to model dependency among neighbouring sites using the stochastic partial differential equation (SPDE) approach (Lindgren & Rue, 2015).Using a penalised complexity prior for the SPDE model (Fuglstad et al., 2019), we set the prior probability that the median range was greater than 50 km 2 to 0.5 and the prior probability that the marginal standard deviation was greater than one to 0.5 (Text SI.1).The spatial relationship between sites was defined by creating a weights matrix, derived from a mesh bounded by the coastline of GB, and the unique locations of all surveys conducted between 1998 and 2016.The model was specified as above (Eq.(2)), with the addition of the spatially structured random effect at survey i, denoted by S i : We then specified a third INLA model, including environmental variables known to influence occupancy; referred to as the "Base model + SPDE + environment".Per species, the impact of including each environmental covariate on the model was assessed using Watanabe-Akaike information criterion (WAIC).We also explored fitting the models with interaction terms between landcovers and the year effects but found no improvements to the models, according to WAIC.Ultimately, model selection resulted in a different set of environmental covariates per species (Table 1).All variables were scaled around their respective means before inclusion in the models (mean = 0, sd = 1).Environmental covariates were included with the survey specific covariates in the matrix X i (Eq.( 3)).
We report the estimated trends as an index of relative occupancy, like the current format of the reported GAM trends (Bat Conservation Trust, 2021).For all models per species, posterior estimates for each year were rescaled so that 1999 was the baseline year, assuming 100% relative abundance.The posterior estimates were sampled 10,000 times to calculate to posterior predictive distribution.The 95% credible intervals around the mean posterior estimates per year were calculated from the posterior predictive distribution and were rescaled so that the credible intervals in the baseline year were zero.

Model comparisons
To quantify the impact of including the SPDE and the environmental covariates on model performance per species, WAIC and model mean absolute error (MAE) were used.To assess the impact of including the SPDE and environment on the change in the index of relative occupancy between the baseline and final year of the time series, the distribution of the posterior estimates for the final year of the time series (2016) were compared between models.We also conducted k-fold cross validation to assess the fit of each model to the data.The data were split based on 11 British regions (Table SI.4), or 10 regions for the E. serotinus models.When implementing INLA using the R-INLA package predictions on missing values were automatic, and per model, per species, data from each of the 11 or 10 regions were replaced with missing values.We evaluated the impact of removing each region on the estimated population trends.

Comparing disaggregated trends
We separately fitted models for Scotland (number of surveys = 565) and England (number of surveys = 5141) for P. pipistrellus and P. pygmaeus.Due to the small number of surveys conducted in Wales (number of surveys = 226), it was not possible to construct stable models for this country alone.We excluded N. noctula and E. serotinus from the Scottish analysis, as the data were highly zero inflated for these species and they are rare or absent in Scotland.We used the best model per species according to WAIC and MAE, and the certainty in the posterior estimates.The posterior distributions of the disaggregated trends were compared to the GB trends.

Impact of spatial and environmental bias on relative occupancy
Trends in relative occupancy indicate in 2016 the populations of P. pipistrellus, P. pygmaeus and N. noctula increased in GB since 1999, whereas there was no change in the population of E. serotinus.Accounting for spatial bias in sampling resulted in small changes to the Pipistrellus spp.trends and improved the model fit for all but E. serotinus.The same pattern was found when both spatial and environmental biases were included compared to the Base model.Under the Base model, the increase in the index of relative occupancy was greatest for P. pipistrellus (160.6%,95% CI 151.1-171.0%).The increases in the index of relative occupancy of the other three species were lower, with larger 95% credible intervals P. pygmaeus (117.6%,95% CI 108.3-128.4%),N. noctula (11.9%, 95% CI 12.5-123.1%%),and E. serotinus (102.96%,95% CI 88.46-120.2%)(Fig. 2; Table SI.5-8).
Accounting for spatial autocorrelation in the data, using the Base model + SDPE resulted in marginal changes to the trends of each species when compared with the Base model trends.The shapes of the trends were broadly the same (Fig. 2).The magnitude of increase in relative occupancy between 1999 and 2016 was slightly greater for all species when estimated using the Base model + SPDE.For instance, the populations of P. pipistrellus and P. pygmaeus increased by 163.0%(95% CI 152.2.6-174.4%)and 18.9% (195% CI 109.0-130.3%),respectively.The joint posterior estimates for 2016 were largely overlapping for both Pipistrellus spp.indicating the estimates from the Base models and Base models + SPDE were virtually identical; there was an 0.58 probability that the estimate for P. pipistrellus was greater using the Base model + SPDE compared to the Base model, and 0.55 for P. pygmaeus (Table SI.9).Over time the estimated trends were also generally higher than those produced using the Base model (Fig. 2; Table SI.5-8).Including environmental covariates (i.e., Base model + SPDE + environment) resulted in similar estimated changes in relative occupancy between 1999 and 2016 to the Base model for P. pipistrellus (159.8%,95% CI 148.2 -172.2%) and P. pygmaeus (116.3%, 95% CI 106.6 -128.0%)(Table SI.5  & 6).The joint posterior estimates per species for 2016 using the Base model + SPDE + environment were largely identical to the Base model and Base model + SPDE estimates (Fig. 2; Table SI.9).
Responses to environmental covariates included per species varied, we defined responses as strong where the 95% credible intervals did not cross zero and the posterior probability was within 5% of zero (negative impact) or one (positive impact).Effects were considered moderate if the posterior probability was within 10% of zero or one.P. pipistrellus was strongly negatively impacted by greater mean spring precipitation.P. pygmaeus was strongly positively impacted by greater mean summer temperatures and mean spring precipitation, but strongly negatively impacted by greater mean winter temperatures.The relative occupancy of P. pipistrellus and P. pygmaeus were strongly negatively impacted by higher percentages of grassland, agriculture/cropland and urban landcover.P. pipistrellus was also strongly negatively impacted by higher broadleaved woodland cover and greater needle-leaved woodland cover moderately negatively impacted P. pygmaeus relative occupancy (Fig. 3; Table SI.10 & 11).Higher percentage cover of broadleaved woodland and grassland moderately negatively impacted the relative occupancy of N. noctula (Fig. 3; Table SI.10 & 11).

Model comparisons
According to WAIC, the models for all species, except E. serotinus, were improved by including the SPDE and environmental covariates (Table 2  Disaggregating trends to country level indicated that the GB trends were primarily representative of the English trend, as is also seen for the country-level GAM trends published by the NBMP (Bat Conservation Trust, 2021).The index of relative occupancy since 1999 in 2016 in England indicated increases for P. pipistrellus (160.1%,95% CI 150.1-170.54%)and P. pygmaeus (+107.5%,95% CI 96.2-119.7%)respectively (Fig. 4, Table SI.13),similar to the British estimates (Fig. 2).Although, the 95% credible intervals of P. pygmaeus for 2016 crossed 100 and the probability that the index of relative occupancy was greater in 2016 in GB than in England was 0.78 (Table SI.14).Trends in Scotland diverged from the British and English trends, with a smaller magnitude of increase in P. pipistrellus (136.6%,95% CI 101.7-178.4%)and a probability that the GB and English indexes were greater than the Scottish indexes in 2016 were 0.75 and 0.76, respectively (Table SI.14).In contrast, there was a greater magnitude of increase in P. pygmaeus between 1999 and 2016 (171.1%,95% CI 131.47-218.7%;Fig. 4; Table SI .13 & 14).This reflected the Scottish NBMP GAM trends (Bat Conservation Trust, 2021), although the population trend for P. pygmaeus in Scotland showed greater fluctuation over time.The 95% credible intervals were wider for the Scottish estimated trends compared to the English and British trends, indicating less certainty (Fig. 4 & Table SI.13).N. noctula and E. serotinus trends in relative occupancy for England broadly followed the British or England-Wales trends.The magnitude of change over time was slightly lower for the English N. noctula trend (105.57%,95% CI 96.75-115.77%)and the English E. serotinus trend was the same as the England-Wales trend (Fig. SI.4;Table SI.14).
Broadly, the direction of impact of environmental covariates remained the same in disaggregated models as for the British wide models (Fig. 4; Fig. SI.4).The 95% credible intervals were generally larger for the Scottish P. pipistrellus and P. pygmaeus models (Fig. 4) and some differences in responses to environmental covariates were found.The relative occupancy of P. pipistrellus was strongly negatively impacted by a higher percentage of urban land in England, but strongly positively impacted in Scotland.P. pygmaeus relative occupancy was moderately negatively impacted by a higher percentage of needleleaved woodland or grassland in England but no there was no impact in Scotland.In contrast, broadleaved woodland strongly positively impacted P. pygmaeus in Scotland but no impact was found in England.Climatic variables had no impact on either Pipistrellus spp. in Scotland (Fig. 4 & Table SI.15).

Discussion
Spatial biases are a common problem in many biodiversity datasets (Buckland & Johnston, 2017;Isaac & Pocock, 2015), particularly those collected by citizen scientists where participation is typically greater in densely human populated areas or where access to sites is easier.Accounting for these biases when modelling population trends is important to ensure trends are robust and representative, especially when trends are used to inform policy decisions at smaller scales than those of the reported trends.Here, we showed that whilst the population trends of four bat species, P. pipistrellus, P. pygmaeus, N. noctula and E. serotinus, were broadly robust to spatial autocorrelation in the data, model performance was improved by including a SPDE for all species, except E. serotinus.There were small changes in the magnitude of change over time when biases in sampling and environmental correlates of occupancy were included in the models.However, these small differences could propagate over time, resulting in very different estimates and ultimately impacting policy decision making.
Nationally disaggregated trends are fitted for England and Scotland where possible from the NBMP Field Survey data by the Bat Conservation Trust (Bat Conservation Trust, 2021).We found spatial biases due to sampling effort still exist in nationally disaggregated data and need to be properly accounted for.We found differences in nationally disaggregated trends compared to the British trends, which is consistent with previous modelling of these data (Bat Conservation Trust, 2021).English trends mostly followed the British trends, whereas Scottish trends showed differences in patterns over time.The magnitude of increase trend was greater for P. pygmaeus in Scotland, indicating this species may be doing better here than in the rest of GB, whereas the change in P. pipistrellus occupancy was lower.These results indicate the British trend masks substantial variation among countries, even when spatial and environmental biases in survey effort are accounted for.It was not possible to construct trends for Wales for any species due to a paucity of surveys in this country.This highlights the poor representation this country has in the British overall trends and action to increase the number of surveys here is required.
The inclusion of environmental variables generally improved the model fit, although did not impact the estimated trends in relative occupancy.Broadly, we found responses to environmental variables supported existing knowledge of habitat and climatic correlates of occupancy for the four species within GB (Andrews et al., 2016;Davidson-Watts et al., 2006;Dietz et al., 2009;Mackie & Racey, 2007).P. pipistrellus and P. pygmaeus relative occupancies were varying impacted by climatic conditions.Positive responses to higher annual mean summer temperatures were identified in England.Contrasting responses to high annual mean spring precipitation interspecifically and internationally and negative responses to higher winter temperatures suggesting projected future climate changes (Betts & Brown, 2021) may have complex effects on these species.Key risks are phenological mismatches between bats and their prey or reduced prey availability due to climate change (Andrews et al., 2016;Froidevaux et al., 2017).Negative responses to a higher proportion of urban cover at the GB level by P. pipistrellus and P. pygmaeus indicate these species are not as urban tolerant as reported in some of the literature (Gili et al., 2020;Jung & Threlfall, 2018).Anticipated urban expansion within GB may cause the current positive trends in these species to be reversed (Sachanowicz  et al., 2019).
Differing responses to environmental variables between England and Scotland emphasised that country level analysis is important in understanding impacts of environment on relative bat occupancy.However, beta coefficient credible intervals were generally large for the Scottish models due to the smaller sample size, which meant strong responses could not be identified.Including an interaction between country as a factor and the environmental covariates may enable country-level responses to be identified.There were unexpected responses to some environmental variables and we suggest treating the responses with caution.For instance, the neutral or negative responses of P. pipistrellus, P. pygmaeus and N. noctula to broadleaved woodland, used by these species for foraging or roosting (Boughey et al., 2011;Mackie & Racey, 2007), may due to the finer scale habitat detail, such as trees within hedgerows, being missed at the 300 m 2 scale of the original landcover data (ESA, 2017).Improved understanding of the impact of local environment on bat occurrence may be achieved by using the exact locations of points or walks along the 3 km transect and finer scale habitat information.

Conclusions
Importantly, we show that spatial biases were impacting the estimated population trends of the species monitored here, although the impact was small.The NBMP Field Survey is a well-structured citizen science survey, but many others rely on more opportunistic data collection or volunteer chosen sites.Spatial biases are a common problem in citizen science datasets (Buckland & Johnston, 2017;Isaac et al., 2016), and the magnitude of impact on estimated population trends for many monitored groups is likely to be greater than the small effect sizes reported here.The models presented here would provide a robust way to ensure spatial biases are accounted for, ensuring representative population trend estimates.Including environmental correlates of occupancy in the models was also important in improving model fit for three species.Sizeable gaps in survey effort remain in areas of lower human population density, such as Wales and the north of England.Finally, whilst the positive or stable trends may be cause for hope, the baseline for these trends is relatively recent.The true extent of historic declines in bats in GB are unknown, but are thought to have been significant, and the populations of these species may still be far below their historic levels.Determining the extent of past declines would be an important step in understanding the significance of these trends.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Locations and number of National Bat Monitoring Programme Field Surveys carried out per country per year.A map of Great Britain shows the locations of the NBMP Field survey sites (n = 552), indicated by black dots (A).The histogram shows the total number of NBMP Field Surveys per country per year (Great Britain n = 5932; England n = 5141; Scotland n = 565; Wales n = 226) (B).

Fig. 2 .
Fig. 2. Trends in relative occupancy over time of P. pipistrellus, P. pygmaeus, N. noctula and E. serotinus.The change in an index of relative occupancy is from a baseline (1999 = 100; horizontal dotted line) produced using presence/absence at 12 points/walks per survey for A) P. pipistrellus, B) P. pygmaeus, C) N. noctula and C) E. serotinus.Three modelled trends are shown per species, estimated using a Base model (black line), including the stochastic partial differential equation (SPDE) (red line) and including the SPDE and environmental variables (blue line).The solid lines show the mean estimated trend and the dashed lines the upper and lower 95% credible intervals around the estimates for each year.

Fig. 3 .
Fig. 3. Effects of site level land-cover and climatic conditions on the relative occupancy of P. pipistrellus, P. pygmaeus, N. noctula and E. serotinus.Points show the mean posterior beta coefficients of environmental covariate fixed effects for models of relative occupancy for A) P. pipistrellus, B) P. pygmaeus, C) N. noctula and D) E. serotinus.Bars show the 95% credible intervals, where these do not cross zero (vertical dotted line) the effects are considered strong and indicated by asterisks.

Fig. 4 .
Fig. 4. Nationally disaggregated population trends in relative occupancy and covariate impacts on occupancy.Comparison of trends in relative occupancy (A and B) and covariate beta coefficients (C and D) fit for Great Britain (orange) and disaggregated to country level, England (red) or Scotland (blue) for P. pipistrellus (A and C), and P. pygmaeus (B and D).Trends show a percentage change in an index of relative occupancy from the baseline year (1999 = 100; horizontal black dotted line).Solid lines show the mean estimated trend and the dashed lines the upper and lower 95% credible intervals.Covariate beta coefficient plots show the mean estimate (dot), 95 % credible intervals (vertical lines), where these do not cross zero (horizontal dotted line) the effects are considered significant.

Table 1
Environmental variables included per species in the Base model + SPDE + environment.The percentage coverage or mean and standard deviation (sd) of each variable across Great Britain (GB) and at the National Bat Monitoring Programme Field Survey (NBMP FS) sites are shown.
E.Browning et al.
; Fig. SI.2).The most competitive model was the Base model + SPDE + environment for P. pipistrellus and P. pygmaeus.The Base model + SPDE and Base model + SPDE + environment were comparably competitive for N. noctula, whereas the Base model was the most competitive for E. serotinus (Table 2; Fig. SI.2).The MAEs were marginally lower for the Base model + SPDE + environment for the Pipistrellus spp.and for N. noctula it was equal to the Base model + SPDE.The MAE of the Base model was lower than the models including the SPDE and environmental covariates for E. serotinus (TableSI.12).K-fold cross validation across models indicated no model was over-fitting to the data.For all species the population trends varied slightly 3.3 Comparing disaggregated trends

Table 2
WAIC values for each species' models.Best performing models where the WAIC was lowest are indicated by an asterisk (*).Models within five WAIC of each other were considered competitive.