Assessing factors associated with changes in the numbers of birds visiting gardens in winter: Are predators partly to blame?

Abstract The factors governing the recent declines observed in many songbirds have received much research interest, in particular whether increases of avian predators have had a negative effect on any of their prey species. In addition, further discussion has centered on whether or not the choice of model formulation has an effect on model inference. The study goal was to evaluate changes in the number of 10 songbird species in relation to a suite of environmental covariates, testing for any evidence in support of a predator effect using multiple model formulations to check for consistency in the results. We compare two different approaches to the analysis of long‐term garden bird monitoring data. The first approach models change in the prey species between 1970 and 2005 as a function of environmental covariates, including the abundance of an avian predator, while the second uses a change–change approach. Significant negative relationships were found between Eurasian Sparrowhawk Accipiter nisus and three of the 10 species analyzed, namely house Sparrow Passer domesticus, starling Sturnus vulgaris, and blue tit Cyanistes caeruleus. The results were consistent under both modeling approaches. It is not clear if this is a direct negative impact on the overall populations of these species or a behavioral response of the prey species to avoid feeding stations frequented by Sparrowhawks (which may in turn have population consequences, by reducing available resources). The species showing evidence of negative effects of Sparrowhawks were three of the four species most at risk to Sparrowhawk predation according to their prevalence in the predator's diet. The associations could be causal in nature, although in practical terms the reduction in the rate of change in numbers visiting gardens accredited to Sparrowhawks is relatively small, and so unlikely to be the main driver of observed population declines.

& Rickard, 2002). Additional evidence in support of this comes from landscape-scale experimental work, delivering supplementary seed to winter farmland, which produced a positive effect on House Sparrow population trends on the study sites (Siriwardena et al., 2007).
The factors behind the declines seen in urban House Sparrow populations, for example, have proved far more difficult to identify and resolve, prompting significant debate. Lack of suitable nest sites (Chamberlain, Toms, Cleary-McHarg, & Banks, 2007;Shaw, Chamberlain, & Evans, 2008), loss of invertebrate food supplies (Peach, Vincent, Fowler, & Grice, 2008), pollution (Summers-Smith, 2007), disease, and increased predation by cats and Eurasian Sparrowhawk (hereafter Sparrowhawk) Accipiter nisus (Bell, Baker, Parkes, Brooke, & Chamberlain, 2010;Churcher & Lawton, 1987) have all been put forward as potential causal factors. Work by Vincent (2005), Peach et al. (2008), Peach, Mallord, Orsman, Ockendon, and Haines (2013), Peach, Sheehan, and Kirby (2014) and Peach, Mallord, Ockendon, Orsman, and Haines (2015) suggests that food availability may play a role, reducing breeding productivity and lowering postfledging survival rates. Shaw et al. (2008), working in the southwest of England, found that House Sparrows were more likely to be lost from the more affluent parts of cities, suggesting that increased demands for off-road parking and the more managed approach to gardening reduced the habitat available to House Sparrows and, by doing so, altered foraging opportunities and predation risk. The loss of large urban gardens through development and "infilling" was found to have a negative impact on urban House Sparrow abundance (Chamberlain et al., 2007).
One area where there has been particularly vigorous debate has been around the possible role of predation in urban House Sparrow declines. For example, data collected through the BTO Common Birds Census (CBC) and the BTO/JNCC/RSPB Breeding Bird Survey (BBS) have charted a 70% decline in the House Sparrow breeding population within England since 1977 and the species has been placed on the Red List of Birds of Conservation Concern (Eaton et al., 2015;Robinson et al., 2015). Much of this debate has centered on the recovering breeding population of Sparrowhawk, a specialist predator of small birds whose English breeding population has increased by 115% since 1975 (Robinson et al., 2015).
This increase in abundance, brought about through improved breeding success following a decline in the use of organochlorine pesticides (Newton & Wyllie, 1992), was accompanied by a recolonization of its former breeding range (Balmer et al., 2013) and an increased use of urban sites (Chamberlain et al., 2005). This increase spans a similar time frame to that of the declines in House Sparrows, prompting some authors to suggest a causal link (Bell et al., 2010).
Previous studies examining the potential impacts of a recovering Sparrowhawk population failed to find any large-scale effect on breeding abundance (Newson, Rexstad, Baillie, Buckland, & Aebischer, 2010;Thomson, Green, Gregory, & Baillie, 1998). The effect of predators on their passerine prey has generally been assumed to compensate for birds that would otherwise succumb to mortality through other means, the predation being compensatory rather than additive (Newton, 1998). The lack of evidence for an effect on breeding numbers is, therefore, perhaps not surprising. Evidence for a compensatory effect may be found through the analysis of postbreeding numbers and analysis of data collected over the winter may give alternative insight into the scale of predation effects. Perrins and Geer (1980) and Newton, Dale, and Rothery (1997) studied the effect of an increase in a Sparrowhawk population on nonbreeding tits (Paridae) and other woodland species and found the seasonal pattern of mortality was altered, as was the peak in numbers. While these studies highlight local impact, they fail to provide evidence of wider scale patterns, something that can only come from a much larger study.
Isolating the impact of predation on prey populations can be challenging, and it is important to note that the method used to model the effect of environmental covariates on changes in populations of birds can also have an effect on the results. Toms (2009), Jones-Todd, Swallow, Illian, andToms (2017), and Bell et al. (2010) all analyzed GBFS data to test for Sparrowhawk effects on garden birds and found contrasting results. Chamberlain et al. (2009) found no significant effect while accounting for temperature and number of feeding units in the model. Bell et al. (2010), however, found significant negative effects of Sparrowhawk on House Sparrows but failed to account for any additional environmental covariates. Jones-Todd et al. (2017) fitted multispecies spatio-temporal models to GBFS data to jointly model changes in the predator and prey species, but once again additional environmental covariates were not modeled explicitly. In addition, the way in which we use the covariates to model change may influence our conclusions. Newson et al. (2010), for example, proposed the use of a change-change model, as opposed to a standard log-linear model frequently used (e.g., Thomson et al., 1998), in the analysis of predation effects on songbirds. Newson et al. (2010) also failed to find any evidence that House Sparrow breeding population declines were linked to an increasing Sparrowhawk population.
The methods used in this paper aim to avoid these two problems in re-assessing the question of Sparrowhawk effects on songbirds by using an array of possible covariates to explain the observed changes in songbird abundance, while using two alternative model formulations to test for consistency in the results.

| Survey methods
The study in this paper uses an extensive volunteer survey conducted annually by the BTO. The BTO GBFS has been monitoring the numbers of birds visiting private garden feeding stations each winter in the UK since 1970 (http://www.bto.org/volun teer-surve ys/gbfs).
For consistency of comparison with other studies, we utilize data from the onset of the study in 1970 up until 2005. GBFS sites are selected to give a representative range of garden types and spatial distribution across the United Kingdom, although they are selected from existing BTO survey volunteers (see e.g., Chamberlain et al., 2009;Chamberlain et al., 2005 for further discussion). It is therefore important to consider any subsequent analysis with the caveat that these gardens belong to observers with an existent interest in birds, and therefore may not be fully representative of the wider population. Participants note down the weekly maximum number of each species seen feeding on the provisioned food across a 26week period spanning the months October to March each year. Any avian predators preying on the birds visiting the feeding station are also noted. There is an average annual turnover rate of 8% for the sites and years considered in this paper, but replacement sites are selected to be similar in local habitat, size, and location to those that have left the scheme. In total over the 36 years, 693 different sites were monitored, with an average of 174 sites per year (five number summary [53, 138.5, 147.5, 226.5, 306]). This was somewhat skewed by the first 3 years, where fewer sites were monitored and after this all but 1 year saw at least 125 sites monitored. In total over the full period, 6,185 individual site-year observations were included in the analysis.
Here, we consider the impact of a single avian predator, namely the Sparrowhawk, on the 10 species of avian prey considered by Chamberlain et al. (2009) (Table 1). Figure 1 shows the variability in the mean number of the 26 weekly observations that are successfully completed each year, across sites monitored in that year. There is clearly a dip in observations during the late 1970s and early 1980s, which has been followed by much more uniform and consistent surveillance since 1990.

| Analytical methods
Due to the fact that multiple observations are conducted within each site every year, it would require very computer intensive methods to analyze all the raw data. Following the approach adopted by Bell et al. (2010), observations across the 26 weeks were therefore averaged to give a mean of weekly maxima per site in each winter that the survey was conducted at that site. Despite the species selected being some of the most commonly recorded under the GBFS scheme, any given species is absent for many site-by-year combinations. This is particularly the case for more localized or habitat-specific species such as Coal Tit, where up to 20% of the annual average site counts were exactly zero. The data therefore relate to an assumed continuous distribution, bounded below by zero, but also with a discrete mass at zero. Due to some sites being particularly large and/or well provisioned, and thus capable of sustaining large numbers of birds, the distributions are often heavily right-skewed. Careful consideration must therefore be given to the specification of the distributional form of the fitted model.
The Tweedie distributions (Dunn & Smyth, 2005;Jørgensen, 1987Jørgensen, , 1997) are a class of distributions belonging to the exponential dispersion models and can be specified through their mean-variance relationship. A Tweedie distribution with mean μ has variance p , where φ is a positive dispersion parameter and TA B L E 1 Proportion of sites where each species is observed at least once during the survey period and percentage of all conducted weekly counts where each species is observed  1970 1975 1980 1985 1990 1995 2000 2005  12 14 16 18 20 22 24 Year Mean no. of counts p ∉ (0, 1) a real-valued index parameter. This class of distributions is incredibly flexible, shown particularly by the fact that they contain many standard distributions as special cases, such as normal (p = 0), Poisson (p = 2), gamma (p = 3), and inverse Gaussian (p = 4).
For values of p ∈ (1,2), the Tweedie distributions are non-negative continuous with a discrete mass at zero. A greater discussion of these distributions is provided in , but, in essence, the distribution corresponds to a zero-inflated distribution when p ∈ (1,2). The Tweedie distributions offer a unified alternative to other zero-inflated approaches, where the proportion of zeros is commonly modeled through a Bernoulli variable, with a separate distribution for the positive observations. The power mean-variance relationship of the Tweedie distributions is also one that is common to ecological applications (Taylor, 1961). These distributions have been previously used in the environmental sciences to analyze fisheries biomass data (e.g., Candy, 2004;Foster & Bravington, 2013) and rainfall data (e.g., Dunn, 2004;Hasan & Dunn, 2010) but not for averaged data such as those analyzed in this paper.
The modeling process adopted in this paper approaches the question of Sparrowhawk effects on each of the 10 prey species independently, using two different methods. First, we use the methodology described in  and extended in , modeling change in the log expected count as a function of spatially varying environmental covariates, some of which also vary over time. Let y it be the mean of weekly maxima at site i in year t for a given species of prey, then the model is expressed in the following terms: where x i is a vector of site-dependent covariates and v it are both siteand time-dependent covariates with regression parameter vectors β and γ, respectively. We refer to this as the "standard model." In order to make full use of the data available, μ i0 (i.e., the expected value at site i in the year before the survey commences at that site), is estimated using a data augmentation approach.
We use the following covariates that depend on site only: northing and easting values obtained from a six-figure grid reference for each site and a two level factor variable denoting whether the site is rural (−1) or suburban/urban (+1). We also include time-varying covariates, namely a year-lagged expected count of the prey species being modeled to test for density dependence, the mean annual count of Collared Dove (except in the analysis of Collared Dove) and Sparrowhawk and number of ground frost days over the months the survey is conducted (October-March) from the Met Office UKCP09 gridded datasets (http://ukcli matep rojec tions.metof fice.gov.uk/). Collared Dove is included as a "pseudo-predator," which has previously been used to test for spurious correlations that may arise coincidentally (Newson et al., 2010;Thomson et al., 1998). The method used to test for density dependence assumes a first-order Markov property that the number of birds observed in year t depends only on the previous year, and is restricted to be negative (Dennis & Taper, 1994).
The ε i denote site-specific random effects such that which account for variation that is specific to the sites but unexplained by the fixed effects. Using hierarchical centering, we can reparameterize Equations 1 and 2 as follows: and We use this parameterization in order to reduce correlation between parameters and improve the efficiency of the parameter estimation process.
The second modeling approach, the "change model," modifies the above methodology into the framework outlined by are used in both modeling approaches and are specified in Table 2.
Covariates are normalized to ensure all covariates are on the same scale and hence priors on the corresponding coefficients are also specified on a meaningful scale (King, Morgan, Gimenez, & Brooks, 2010). We are particularly interested in which environmental factors best explain the observed changes in populations of songbirds.
We therefore include an array of environmental covariates and use Goodness-of-fit is assessed through the Bayesian p-value p , a measure of posterior predictive fit (Gelman et al., 1996). The deviance is used as the discrepancy statistic and p-values outside the interval 0.025 ≤ p ≤ 0.975 would give rise to evidence of poor fit of the model.

| RE SULTS
The two models were run independently for all 10 species. The standard model converged very quickly. Using 20,000 iterations with the first 5,000 being discarded as burn-in appeared to be very conservative. The change model took longer to converge and here we ran 100,000 iterations with the first 60,000 iterations discarded as burn-in. Posterior summary statistics from the two analyses are presented in Tables 3 and 4. For the regression parameters, significance was assumed for covariates with Bayes factors >3. Goodness-of-fit was assessed through Bayesian p-values using the deviance as the discrepancy statistic (Table 5). The only analysis giving possible concern was Coal Tit under the standard model, with an estimated pvalue just inside the rejection region. This may be due to a marginally bimodal distribution for this species which the Tweedie distributions have difficulty in fitting to.

| Sparrowhawk effect
Across the two models, eight of the possible 20 correlations between prey species and Sparrowhawks were found to be significant.
Of these, five came from the standard model and three from the Gotmark and Post (1996) estimated the relative predation risk for a number of species of songbirds to predation by Sparrowhawks, a measure of how frequently the species occurs in the Sparrowhawk's diet having normalized for its prevalence in the surrounding environment. In particular, species frequently found feeding on the ground have been shown to be particularly susceptible to Sparrowhawk predation (Chamberlain et al., 2009;Götmark & Andersson, 2005).

Associations between high relative predation risk and negative
Sparrowhawk coefficients may give further support to a genuine Sparrowhawk impact on prey populations.
Relative predation risks for all of the species analyzed in this paper, aside from Collared Dove, were calculated by Gotmark and Post (1996)  This suggests that those species appearing most at risk to predation by Sparrowhawks are also showing negative relationships with Sparrowhawk abundance. Of the four species most at risk from predation by Sparrowhawk, three had significant negative estimates under both our models. Under the relative predation risk, a value of zero is equivalent to a Sparrowhawk randomly searching for prey.
According to the linear model fitted to our results and predation risk, a zero value of relative predation risk relates to a value close to, but slightly below, zero under the two models used in this paper. This could relate to the fact that a zero relative predation risk does not

| Effects of other variables
With reference to the other regression covariates, easting was found to be a more important predictor of changes in songbird numbers than was northing. Northing was only significant for one species, positive for House Sparrow in the change model. However, Coal Tit and Greenfinch were the only species with significant estimates to show consistency across the two models with respect to these two covariates, both having negative estimates for easting. The other time-invariant covariate, the rural or suburban factor variable, showed some differences between the two models. Where significant, the parameter was always negative suggesting that urban sites are on average faring worse than their rural counterparts.
In the standard model, evidence of density dependence was found in all but two of the species, namely House Sparrow and Estimates of the random effect variance σ 2 were fairly similar across model specifications, suggesting the fixed effects were capturing a similar amount of the variation in both cases. There was, however, a ten-fold difference between the species with the smallest (Robin) and largest (Collared Dove, House Sparrow, Greenfinch) variances. Estimates of the Tweedie parameters φ and p were also consistent under both model formulations. This analysis was conducted on data relating to garden birds attracted to feeding stations. The nature of the survey protocol, as mentioned above, requires that we must be careful in extending inference from these results to wider populations of these species. Chamberlain et al. (2005) found c.96% correlation between winter abundance at GBFS sites and annual indices of relative breeding population change derived from CBC data for both House Sparrow and Starling. Due to the strong correlation in these cases, it seems reasonable to expect that the studied gardens are a good representation of wider populations. Observations of Blue Tits from GBFS were weakly negatively (−0.32) correlated with equivalent breeding density. Nonetheless, the species showing significant negative relationships with Sparrowhawks are also largely those most at risk from predation from Sparrowhawks, appearing more frequently in the Sparrowhawk's diet than would be expected based on their overall abundance (Figure 1).

The significant positive relationships of Sparrowhawk with Coal
Tit and Blackbird in the standard model may well represent confounding factors that lead to Sparrowhawks recolonizing sites that were also more attractive to the former species, rather than a causal relationship. We cannot be certain that the same interpretation does not hold in reverse for the species with negative effects; however, the lack of positive effects under the change model framework and the relationship with predation risk add to the evidence for the negative effects. Sparrowhawks have been encountered previously, may lead to sublethal predation-risk impacts, such as limited access to food and result in increased levels of starvation and a decline in population size (Cresswell, 2008;Seress, Bókony, Heszberger, & Liker, 2011). The reality may be some combination of these two explanations.
The prevalence of significant negative easting effects over northing is most likely due to the higher level of intensive agriculture seen in the eastern UK. This intensification has been shown to be closely linked to negative changes in the populations of many farmland birds (Fuller et al., 1995;Newton, 2004). This might add support to the hypothesis that the intensification of farming has led to decreases in House Sparrow populations. The decreases in farmland may be having an additional effect on the surrounding gardens in this region, something perhaps evident in the regional patterns found in House Sparrow productivity, as revealed by BTO Garden BirdWatch, another study focused on urbanized landscapes (Morrison, Robinson, Leech, Dadam, & Toms, 2014). Coal Tit is a species of coniferous woodland, that makes greater use of garden feeding stations when access to cone crops is limited by a poor seed year (Mckenzie, Petty, Toms, & Furness, 2007). The establishment and maturation of new conifer plantations may influence Coal Tit distribution and we would therefore expect there to be some pattern to the species distribution linked to easting.
The parameter relating to the categorical variable of suburban or rural habitat is, where significant, negative, suggesting that urban populations are faring less well than their rural counterparts. The negative effect of urbanization suggests that there are additional factors common to urban gardens above the level of individual site variation that cannot be explained through the other fixed effects.
Previous studies on rural and urban populations confirm that different trends are taking place in the two different habitats, with populations in urban environments generally doing worse (e.g., Beckerman, Boots, & Gaston, 2007;Newson, Ockendon, Joys, Noble, & Baillie, 2009;Robinson, Siriwardena, & Crick, 2005;Vincent, 2005). The factor variable is only able to pick up general differences between the two habitat types that cannot be explained by the other covariates or the random site effects. Reasons for the significant differences in these subpopulations may warrant further study.
Our results also suggest that density dependence is an important factor governing the numbers of most species observed. It is perhaps not surprising that this is the case given that our data come from garden feeding stations, where birds may well be competing for food within a small area and where other factors linked to the densities of birds present, such as disease transmission (Lawson et al., 2012), may have a role to play.
As outlined in the methods, we included Collared Dove as a "pseudo-predator" to act as a control alongside our investigations  Chamberlain et al. (2005). The use of garden feeding stations by Coal Tits has been shown to be influenced by the size of conifer seed crops, the birds switching to feed on supplementary food more often in years with few cones than in mast years (Mckenzie et al., 2007). Access to the conifer seeds may also be influenced by win- lack of food. The larger effect on counts of birds in relation to these fluctuations in temperature is therefore not too surprising.
We found a fairly large degree of variation in the consistency of species trends across sites. Robins and Blackbirds, for example, seemed to show relatively small intersite differences, while House Sparrow, Starling, Collared Dove and Greenfinch showed less consistency. This suggests that there may be additional factors affecting the latter species that we have not considered in our analyses.
The species with larger random effect variances are ones that tend to form larger groups or flock together. In addition to the factors modeled, there are clearly additional factors that encourage these species to form larger groups at some feeding stations than others, such as the size of site that is not implicitly modeled here.
We see the Tweedie distributions as a widely-applicable but underused tool in ecology and suggest they could become much more widely used outside of their current limited applications. The meanvariance relationship has been shown to be one common to ecology, while remaining incredibly flexible and avoiding the need for strong assumptions to be made about the distribution a priori. Using the two alternative model formulations in general provided consistent results in this case. Overall neither model formulation appeared to consistently outdo the other in terms of posterior predictive power; most variation in Bayesian p-values was between, rather than within, species. We propose that where possible, the use of multiple model formulations is advantageous, allowing different hypotheses and ecological mechanisms to be tested or highlighted. In addition, we have also shown that how the covariates enter a model is not arbitrary. We recommend that careful thought should always be given to the model specification to ensure it is both statistically and ecologically sensible.

ACK N OWLED G EM ENTS
We would first like to thank all of the BTO Garden Bird Feeding Survey volunteers for many years of effort in recording garden birds. GBFS volunteers are recruited from the weekly BTO Garden BirdWatch survey which, through the generosity of its participants, also provides the financial support to keep GBFS running. Particular thanks are due to Clare Simm, Tim Harrison, Andrew Cannon and the late David Glue, who have looked after GBFS, and to Margaret Askew, Fran Bowman, Donna Hobbs, Carol Povey, Alec Prior, Jacky Prior, Heather Pymar and Nicky Ward who have helped with data curation. BTS was part-funded by EPSRC/NERC grant EP/10009171/1.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
BS, STB, and RK conceived the ideas and designed methodology; BS analyzed the data; BS led the manuscript writing with input from MPT. All authors contributed critically to the drafts and gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
The cleaned data file used in this analysis has been placed on Dryad with DOI https ://doi.org/10.5061/dryad.v8j1144.