Reassessing the observational evidence for nitrogen deposition impacts in acid grassland: spatial Bayesian linear models indicate small and ambiguous effects on species richness

Introduction

Nitrogen deposition (Ndep) is a significant threat to the plant diversity of various habitat types, both in north-western Europe (UK National Ecosystem Assessment, 2011) and around the world (Phoenix et al., 2006). The evidence for this position comes from a variety of sources, including correlative analyses of observational data (e.g., Maskell et al., 2010), typically conducted across relatively large areas as “space-for-time” gradient studies, and small scale experiments (e.g., Van der Eerden et al., 1991), although the latter have also frequently been pooled across larger areas through meta-analyses or other approaches to evidence synthesis (Clark et al., 2007; Phoenix et al., 2012; Soons et al., 2017). Reviews of Ndep impacts on plant biodiversity have typically drawn on all of this evidence (Bobbink et al., 2010; UK National Ecosystem Assessment, 2011; Stevens et al., 2011c; RoTAP, 2012; Rowe et al., 2017), in addition to other types of studies, such as Before-After surveys of historic plots (e.g., Britton et al., 2009). Observational and experimental studies are therefore both generally considered useful ways of understanding pollutant-driven biodiversity change in terrestrial ecosystems.

Different inferential approaches are often considered complementary, with large-scale, observational methods potentially allowing access to “treatment” effects across pre-existing gradients, with levels of replication that would likely be challenging to resource via an experimental route (but see Fraser et al., 2013). One cost of this approach is that the effect of interest is likely to be crossed in various complex ways with numerous other variables, including historic drivers for which data are likely to be inaccessible, leaving one with a large choice of covariates that could potentially be included in a model, including some which will be unknown, or suspected to be of importance but impossible to access. Furthermore, spatially autocorrelated variables that are not captured by the covariates included, or other processes causing spatial structure such as dispersal, may also need to be accounted for to ensure accurate estimation of regression coefficients (Beale et al., 2010; Crase et al., 2014). The ultimate purpose of a statistical model must also be taken into account in making analytical decisions: does one primarily wish to make predictions, or is the focus on unbiased effect estimation to develop causal understanding (Mac Nally, 2000; Stephens, Buskirk & Del Rio, 2007; Shmueli, 2010)? Although causal inference may imply predictive success, models that are constructed using methods that solely seek to maximise predictive accuracy will not necessarily capture causal processes accurately. Even when explicitly aiming for causal explanations through regression modelling in non-experimental situations, estimated coefficients may still only have a weak claim to be viewed as causal effects (Gelman & Hill, 2007; Young, 2018). Statistical issues then, in addition to domain-specific understanding, must also be at the forefront when attempting to make statements about cause and effect from observational data (Rubin, 2004; Young, 2018).

The work on Ndep presented here arose from a desire to use information from existing studies to inform the analysis of new data through the use of informative priors in a Bayesian framework (Lemoine, 2019). However, causality-focused estimates of regression coefficients from observational studies of Ndep impacts on plant richness proved hard to find in the current literature. For example, several studies have presented “headline” estimates from univariate models after multiple regression modelling (Stevens et al., 2004; Maskell et al., 2010; Field et al., 2014), implying regression coefficient values for Ndep that do not necessarily have any causal meaning, and/or have used P-value-based stepwise variable selection (PSVS) to derive final models, a strategy long known to be poorly suited to the accurate estimation of regression coefficients for inferential purposes (Greenland & Neutra, 1980; Mac Nally, 2000; Whittingham et al., 2006; Faraway, 2014; Harrell, 2015; Heinze, Wallisch & Dunkler, 2018). Re-analysis of existing studies with the primary aim of developing models that focus on causal inference should therefore be valuable in exploring the dependence of earlier conclusions on modelling choices. Fortunately, data focusing on UK acid grasslands from two such studies (Stevens et al., 2004; Maskell et al., 2010) were available for re-analysis; but note also that very similar work has been done across other habitats (e.g., Maskell et al., 2010; Field et al., 2014), and across larger areas (e.g., Dupré et al., 2010).

The main focus of the current work is on deriving estimates of the effect sizes of nitrogen deposition on plant species richness, rather than on their statistical significance (Amrhein, Greenland & McShane, 2019), whilst accounting for previously unmodelled broad-scale spatial autocorrelation. The model forms investigated were specified a priori; in general these were delimited by the full sets of covariates previously investigated by the original studies, given that these all have good ecological reasons for inclusion. We do not make the (unprovable) claim that our “full” models are fully correct with respect to reality, merely that the inclusion of as many plausible “pre-treatment” covariates as possible is likely to help avoid omitted variable bias, given that the ignorability assumption (i.e., all confounders are measured) is more likely to be satisfied (Gelman & Hill, 2007; Young, 2018), and that even our worst case ratio of response data points to predictors meets rules of thumb put forward by statisticians to help ensure that coefficients can be estimated accurately (Harrell, 2015, p. 72; Heinze, Wallisch & Dunkler, 2018). In addition, we also include a model for the Stevens et al. (2004) data that focuses on a reduced set of covariates chosen for their expected similarity to the ecological impacts of the covariates used by Maskell et al. (2010) for comparative purposes. Two additional models for the Maskell et al. (2010) data were specified post hoc, due to a desire to investigate differences between our results and those of the original paper.

Note that some covariates of interest (e.g., the topsoil variables pH, Al, C:N, and %N) could themselves be influenced by Ndep, but will also have independent impacts on the dependent variable (species richness). These types of partly intermediate variables are sometimes distinguished from fully post-treatment variables as “proxy variables” (Angrist & Pischke, 2009, pp. 64–68) whose inclusion is often better than their omission when a causal interpretation is desired. We assume here that, due to their pre-treatment importance for richness, adjusting for these covariates is more likely to result in accurate estimates of the effect of Ndep than not adjusting for them (Angrist & Pischke, 2009); indeed, in some cases, their inclusion will be essential for avoiding bias (Rosenbaum, 1984). Ultimately, we assume that these potential post-treatment (or proxy) variables are “plausible surrogate[s] for … clearly relevant but unobserved pretreatment variable[s]” (Rosenbaum, 1984).

Materials and Methods

Results

The inclusion of both the 1 km square random effect and the spatial random field improved almost all the models of the effects of N deposition on plant species richness in acid grassland explored here (as judged by DIC, see “Statistical models” above). The one model (SEA04 model 1 with bryophyte data) where ΔDIC was negative when either the 1 km square random effect or the spatial random field was dropped (indicating that the simplified model was preferred) still required a spatial random effect at some scale: dropping both random effects gave a DIC value of 1732.7 compared to a value of 1712.7 for the full model; the preferred model in this case was that with the spatial random field only (DIC = 1709.0). Estimates of the overall random effects variance explained by the structured spatial field varied between 87.4% (MEA10, vascular plants only) to 99.7% (SEA04 Model 1, vascular plant and bryophyte richness combined) across the models including both random effects. All estimated spatial random fields can be viewed in Supplementary Information 1.

The addition of an extra kilogram of total N deposition on plant species richness per hectare per year was consistently amongst the smallest effect estimated by our models, with low uncertainty in comparison to other covariates (Figs. 1A–1C; Supplementary Information 2: Figs. A2.1, A2.2, A2.3). Exponentiating the fixed effects indicated that the impact of increasing total N deposition is equivocal across these datasets and models (Table 2): six models estimated a drop of around 1% in species richness as the most likely impact (Figs. 1B–1C; Supplementary Information 2: Figs. A2.2, A2.3; Supplementary Information 3: Figs. A3.1, A3.2), although two of these models (Supplementary Information 3) were formulated post hoc and were not favoured as the best models for the dataset by their DIC estimates. Two models, estimated from MEA10 using the DIC-favoured model, suggested that an increase in species richness had more support (Fig. 1A; Supplementary Information 2: Fig. A2.1; Table 2). The reduced covariate set used for the SEA04 data (model 1) for closer comparison to the covariate set used in the MEA10 analysis made no difference to the estimated effect of N deposition calculated from the larger set of covariates from SEA04 (model 2; Figs. 1B–1C; Supplementary Information 2: Figs. A2.2, A2.3; Table 2). Figure 2 illustrates the estimated partial effect of N deposition (with other covariates set to zero) on vascular plant species richness for two of our models (MEA10 and SEA04 model 1); these demonstrate both the relatively small predicted mean impact on species richness, the disagreement in impact direction between the datasets, and the uncertainty attached to these predictions.

The results presented in SI3 were included as exploratory, post hoc, analyses after it was found that our DIC-favoured models of the MEA10 dataset indicated that the most likely overall influence of Ndep on species richness was more likely to be positive, contrary to the conclusions of Maskell et al. (2010). Our additional analyses (SI3), each with only one of the spatial random effects used in the full model, did show a reversal of average sign, with the analysis matching that of Maskell et al. (2010), i.e., using a 1 km square random effect, indicating the strongest negative effect of N deposition on vascular plant richness (a mean of a 2% loss of species richness per kg ha⁻¹ yr⁻¹ of total N deposition; Table 2; Supplementary Information 3: Fig. A3.1). However, this model had a DIC value 7.4 units higher than the favoured model.

Discussion

Previous work on the estimation of the effects of nitrogen deposition on plant species richness from space-for-time substitution studies is likely to have had a significant impact on both science and policy (assuming that the number of citations received can be used as an index to this). Therefore being clear about the reliability and accuracy of estimates from observational work is important, not only for fundamental reasons of improving ecological understanding and analysis, but also because policy decisions must be made concerning the utility of funding work in this area relative to other ways of assessing environmental change (although field experimental work may have issues of its own; Peters, 1991; Phoenix et al., 2012). Clear and accurate estimates of the likely sizes of effects are invaluable for informing future studies, for example, through power analyses, or the use of informative priors for Bayesian approaches. Our results suggest that only small effects of nitrogen deposition on species richness may be detectable in these observational datasets.

The headline result of Stevens et al. (2004) reported “a reduction of one species per 4-m² quadrat for every 2.5 kg N ha⁻¹ yr⁻¹ of chronic nitrogen deposition”. An estimated reduction in species richness of 23% (based on a change from deposition of 5 kg N ha⁻¹ yr⁻¹ to 17 kg N ha⁻¹ yr⁻¹) was also highlighted as a key message using a regression coefficient estimated from a model including only Ndep as a covariate. The estimate from our SEA04 model 1 implied a loss of around 1% of richness per kg N ha⁻¹ yr⁻¹; however, this ignores the uncertainty associated with this particular coefficient estimate, which included an estimate of no effect at the 97.5% percentile, and the larger uncertainty associated with alternative model specifications (including ones that we have not explored here). Indeed, our estimates from the larger, random-stratified national survey of Maskell et al. (2010) suggested that a small average increase in richness along the total nitrogen deposition gradient was more plausible under the favoured model for that dataset. Furthermore, the partial effects calculated from our full models (Fig. 2) indicated that the average effect on richness across Britain may only approach a difference of around one species in either direction. This is also notable because previous commentaries (e.g., Tipping et al., 2013) have suggested that differences in the apparent effect sizes of Ndep between the two observational studies reanalysed here were likely due to study design, whereas our work suggests that the differences are actually minimal, and that those found previously may simply be artefacts of the different statistical modelling procedures adopted by the original studies. The different median signs found in our analyses could also be due to differences in the statistical populations targeted, namely the differential sampling of grassland types, with the greater plant community coverage of the dataset of Maskell et al. (2010) potentially covering more sites on the left of a net primary productivity/diversity unimodal curve, resulting in an average increase in richness under fertilisation (see below for more discussion of the potential impacts of non-linear responses and other issues on the current work).

Our results should cause others to re-evaluate their approaches to observational data. For example, Field et al. (2014) also used PSVS, and other P-value based selection techniques, to formulate models in their space-for-time analyses of Ndep impacts for a number of semi-natural habitats across Britain, as well as highlighting simple univariate relationships through scatter plots. These authors, however, did note in their methods section that their results should be “interpreted with caution”, although this statement of uncertainty was not clearly carried through to other parts of their paper, nor to other research. For example, Payne et al. (2017) used data from both SEA04 and Field et al. (2014) to forecast the impacts of Ndep on plant species richness under future nitrogen deposition scenarios. Ndep and its polynomial transformations were the only covariates in these models. Payne et al. (2017) noted in their online supplementary material (their Web Panel S1) that such forecasting makes the strong assumption that covariances between variables remain the same at the future time point (i.e., the distribution of nitrogen deposition amounts will vary with, e.g., mean maximum July temperature in the same way in 2030 as it did when originally modelled). However, they did not clearly acknowledge that the Ndep-only models in their forecasting exercise had received no validation in support of their status as the best predictive (or explanatory, for that matter) models by the original studies or elsewhere. Ndep may be a useful covariate for producing good predictive models of species richness in certain habitats, but this should be demonstrated using measures of out-of-sample predictive ability, and the resulting partial effects of Ndep in such models have no intrinsic claim to be reliable estimates of its causal relationship with richness (cf. Payne et al., 2017).

These issues from the literature should not detract from the fact that there are also several ways in which the work presented here could potentially be improved upon. Our efforts are linear models (as for SEA04, MEA10, Field et al., 2014, among others), but there is evidence that the response of species richness to Ndep may be better modelled as unimodal (Tipping et al., 2013; Simkin et al., 2016; Clark et al., 2019). Simkin et al. (2016) still, however, reported declines in the plant species richness of open habitats above 8.7 kg N ha⁻¹ yr⁻¹, very similar to the 7.9 kg N ha⁻¹ yr⁻¹ richness decline threshold identified for acid grassland by Tipping et al. (2013), suggesting that a large part of the gradient studied by SEA04 and MEA10 may still be well-approximated by a linear relationship (although we note that the approach of Tipping et al. (2013) also makes the assumption that there are no omitted linear or non-linear variables correlated with Ndep that might change their estimated univariate breakpoint relationship). Indeed, the linear component of a non-linear trend is often considered the most policy-relevant summary by those who routinely produce ecological indicators (e.g., Soldaat et al., 2017). Different approaches to accounting for this likely non-linearity in multivariable models, such as the inclusion of interactions or smoothers, could however be further explored, particularly given that novel methods focused on causal inference that can account for these issues continue to be developed (Dorie et al., 2019).

Other uncertainties relating to our conclusions pertain to the fact that the broad-scale spatial field used here may be accounting for information that, if known, would change the size of the Ndep impact regression coefficient. This could be in the form of additional covariates, or more highly resolved estimates of the Ndep load that a location has actually received. The Ndep estimates used here (and by the original studies) are resolved to a grain size of 5 × 5 km, and this additional uncertainty could have attenuated our estimate of the regression coefficient (the measurement error in explanatory variables problem; e.g., see Fox, 2016), even if there is no systematic bias relating to the association of 5 × 5 km-estimated Ndep levels with particular types of vegetation. This argument is likely to apply to many of the variables used here and in other studies, particularly given that many are estimated from other modelling exercises at relatively large spatial grain sizes (Table 1). Conversely, measurement error coupled with conditioning on statistical significance, as happens through PSVS (e.g., Stevens et al., 2004; Field et al., 2014), is likely to result in the overestimation of effects (Gelman & Carlin, 2014; Loken & Gelman, 2017). Whether or not the potential for these biases is more serious for inference than the absence of covariates that are unavailable, such as historic land management events (e.g., Rackham, 1986), is difficult to say.

The new observational results presented here are in line with much of the experimental literature on Ndep impacts. For example, Phoenix et al. (2012) reviewed a group of nine experiments conducted across the UK, with N treatments which had been running for lengths of time between 7 and 22 years at the time of review. These studies aimed to examine the impacts of “modest treatment doses and avoid[ed] single dose or solid form applications” (Phoenix et al., 2012) in order to overcome previous criticisms relating to (potentially unrealistic) very high experimental loadings. These nine studies reported no effects of their experimental N treatments on higher plant richness (although these conclusions of no effect were all based on thresholding P-values); two sites indicated increases in richness using the calculated accumulated Ndep dose over the duration of an experiment, although these cases were discounted by Phoenix et al. (2012) as being of either minor ecological significance or transitory. Phoenix et al. (2012) put forward several reasons why the results from this series of experiments may not reflect the true impact of chronic Ndep on plant richness: for example, sensitive species might already have been lost prior to the establishment of an experiment, and/or an experiment might not have been running long enough for the impacts to have been fully realised. Of course, important changes in community structure can also occur without species loss (Hillebrand et al., 2018).

Longer running studies, such as the Park Grass experiment at Rothamsted, England (1856 to the present day; Silvertown et al., 2006), avoid some of these criticisms. Although the Park Grass plots that have received experimental N addition have received very high doses, with annual N fertiliser doses starting at 48 kg N ha⁻¹ yr⁻¹ (Storkey et al., 2015), control plots having received only atmospheric deposition may be a useful comparator for some habitats in the wider landscape. Plot 3 of the Park Grass experiment, for example, has in theory only ever received ambient Ndep (Storkey et al., 2015); Lawes & Gilbert (1880) reported an average of around 48.5 plant species (including bryophytes) in the 0.5 acre (∼2,000 m²) neutral grassland plot between 1862 and 1877; the average between 1939 and 1948 was 34 (Brenchley & Warington, 1958), indicating a decline over this period that preceeded the local increase in Ndep (see the first figure of Storkey et al., 2015 for the local Ndep trend at this location over the 20th century). Unfortunately, the local increase in Ndep at Rothamsted coincided with the decision to split the plots into different liming treatments (Williams, 1974; Storkey et al., 2015) meaning that subsequent richness estimates cannot be unambiguously compared to the earlier numbers (there is confounding of the Ndep increase with the application of lime, and the area that could be directly compared shrinks to a quarter of the original plot, introducing the need for species–area adjustments). Over the period 1991–2012, however, survey data indicate that there may have been a very slight recovery in richness in Plot 3 (although Storkey et al. present no statistics for this trend), coinciding with a reduction in the N composition of the plot herbage and increases in Simpson’s diversity index (Storkey et al., 2015).

Overall, then, experimental data using realistic applications of Ndep appear to support our finding that richness is a relatively insensitive metric of such impacts (see also Hillebrand et al., 2018). Finally, and to avoid any misunderstanding as to the thesis being presented in this paper, we note that we are not stating that overall eutrophication from all sources (e.g., livestock, local fertiliser drift etc.) is unimportant for the conservation of plant biodiversity. General signals of eutrophication in plant communities are widespread and beyond doubt (e.g., Smart et al., 2003b; Willi, Mountford & Sparks, 2005). We are forced to conclude, however, that the contribution of Ndep to this phenomenon appears to be smaller, and more uncertain, than many previous analyses of space-for-time Ndep gradients have concluded (cf. Stevens et al., 2004; Field et al., 2014; Payne et al., 2017). It is possible that this is due to the fact that earlier losses of species due to accumulated historic deposition means that the remaining opportunity to detect effects in recent space-for-time studies is limited; this suggests, however, that richness should not be used as an indicator of Ndep impacts, and previous work taking this approach (particularly where highly significant statistical concerns abound) should no longer be cited in support of general conclusions regarding these impacts.

Conclusions

The evidence for large negative impacts of nitrogen deposition on plant species richness put forward through analyses of observational data appears to have been overstated. We estimate a possible decline in richness of around 1% per kg ha⁻¹ yr⁻¹ of total N deposition from two spatially and temporally separated British space-for-time gradient studies, considerably less than the estimates implied previously by described and/or plotted relationships from primary studies (Stevens et al., 2004; Maskell et al., 2010). Moreover, even this estimate appears uncertain, and our favoured model for the acid grassland data of Maskell et al. (2010) suggests that an average increase in richness of a similar magnitude may be more likely. The previous lack of presented regression coefficients derived with causal inference as their main focus, and of models that account for broad-scale spatial autocorrelation, is important to note, because scientists wishing to use the estimated effects of Ndep for the design of future studies, or for the construction of informative priors in new analyses, may be misled as to the size of effect that is expected to be detectable in datasets of this type. The various models presented here could be thought of as a small section of the “multiverse” of potential approaches to these data, an approach that has been put forward as an additional route to transparency and reproducibility in science, and which can provide insights into the fragility or robustness of particular conclusions (Steegen et al., 2016). As such this work is unlikely to be the last word on these datasets, or in the general area of observational studies of nitrogen deposition impacts. We hope that our re-analyses inspire further efforts to accurately extract the maximum available knowledge from these valuable datasets, whether for explanatory or predictive purposes, and that evidence synthesis in this area takes these uncertainties and methodological issues into greater account going forwards.

Supplemental Information