Advertisement

Hierarchical change-point regression models including random effects to estimate empirical critical loads for nitrogen using Bayesian Regression Models (brms) and JAGS

Open AccessPublished:November 01, 2022DOI:https://doi.org/10.1016/j.mex.2022.101902

      Highlights

      • Hierarchical change-point regression models are suitable for estimating critical empirical loads.
      • The Bayesian framework of these models provides the inclusion of the current critical load and various confounding or modifying variables.
      • Here we present two ways of implementing hierarchical data sets in Bayesian change-point regression models using JAGS and brms.
      • The used model data are foliar N:P ratios of 10 different conifer tree species from 88 European forest sites and 9 different countries covering 22 years (1995-2017).

      Abstract

      The concept of critical loads is used in the framework of the Convention on Long-range Transboundary Air Pollution (UNECE) to define thresholds below which no damaging effects on habitats occur based on the latest scientific knowledge. Change-point regression models applied in a Bayesian framework are useful statistical tools to estimate critical empirical loads. While hierarchical study designs are common in ecological research, previous methods to estimate critical loads using change-point regression did not allow to analyse data collected under such a design. This method update provides an implementation of hierarchical data structure by including random effects such as study sites or as in this example tree species within the Bayesian approach of change-point regression models using two different approaches. The example data set is an European wide gradient study of the impact of climate change and air pollution on forest tree health assessed by foliar nutrient status of nitrogen (N) to phosphorus (P) from 10 different conifer tree species originated from 88 forest sites and 9 countries covering 22 years (1995-2017). Both modelling approaches using JAGS and Bayesian Regression Models using ‘Stan’ (brms) resulted in reasonable and similar estimations of the critical empirical load for nitrogen (CLempN) for temperate forests. These methodological examples of using different approaches of Bayesian change-point regression models dealing with random effects could prove useful to infer CLempN for other ecosystems and long-term data sets.
      • Hierarchical change-point regression models are suitable for estimating critical empirical loads.
      • The Bayesian framework of these models provides the inclusion of the current critical load and various confounding or modifying variables.
      • Here we present two ways of implementing hierarchical data sets in Bayesian change-point regression models using JAGS and brms.

      Graphical abstract

      Keywords

      Method name

      Tabled 1
      Subject areaEnvironmental Science
      More specific subject areaForestry research
      Method nameRevised Bayesian change-point regression models including random effects using brms and JAGS
      Name and reference of original methodChange-point models applied in a Bayesian context by
      • Roth T.
      • Kohli L.
      • Rihm B.
      • Meier R.
      • Achermann B.
      Using change-point models to estimate empirical critical loads for nitrogen in mountain ecosystems.
      Resource availabilityModel data can be found in
      • Du E.
      • van Doorn M.
      • de Vries W.
      Spatially divergent trends of nitrogen versus phosphorus limitation across European forests.
      . Codes including models and graphics are written in R

      R Core Team, Team R Development Core, R: A Language and Environment for Statistical Computing, 2021, http://www.r-project.org/.

      and given as an RMarkdown

      J.J. Allaire, Y. Xie, J. McPherson, J. Luraschi, K. Ushey, A. Atkins, H. Wickham, J. Cheng, W. Chang, R. Iannone, rmarkdown: Dynamic Documents for R, 2021. https://github.com/rstudio/rmarkdown.

      output. The MCMC simulations of the change-point modelswere conducted using JAGS, version 4.3.0
      • Plummer M.
      JAGS: A Program for Analysis of Bayesian Graphical Models using Gibbs Sampling.
      , executed in R using rjags

      M. Plummer, rjags: Bayesian Graphical Models using MCMC, 2019. https://cran.r-project.org/package=rjags.

      and in STAN with brms
      • Bürkner P.-C.
      brms : An R Package for Bayesian Multilevel Models Using Stan.
      .

      Method details

      Bayesian change-point regression model settings

      The revised change-point models using either JAGS or brms are based on the current method of using Bayesian change-point models as in the example given by [
      • Roth T.
      • Sprau P.
      • Naguib M.
      • Amrhein V.
      Sexually selected signaling in birds: A case for Bayesian change-point analysis of behavioral routines.
      ]. The Bayesian analysis is based on MCMC methods [
      • Link W.A.
      • Cam E.
      • Nichols J.D.
      • Cooch E.G.
      Of Bugs and Birds: Markov Chain Monte Carlo for Hierarchical Modeling in Wildlife Research.
      ] with a similar setting used by [
      • Roth T.
      • Sprau P.
      • Naguib M.
      • Amrhein V.
      Sexually selected signaling in birds: A case for Bayesian change-point analysis of behavioral routines.
      ] conducted using JAGS (version 4.3.0 [
      • Plummer M.
      JAGS: A Program for Analysis of Bayesian Graphical Models using Gibbs Sampling.
      ] and executed in R using rjags [

      M. Plummer, rjags: Bayesian Graphical Models using MCMC, 2019. https://cran.r-project.org/package=rjags.

      ]. Posterior distributions were based on parallel chains (t.n.chains=2) with 100000 iterations each (t.n.iter=100000), discarding the first 50000 values (t.n.burnin=50000) and thinning the remainder by 2 (t.n.thin=2). The two parallel chains were used to assess convergence with Gelman and Rubin’s diagnostics [
      • Brooks S.P.
      • Gelman A.
      General Methods for Monitoring Convergence of Iterative Simulations.
      ] using trace and marginal density plots (Fig. 2) with the R function plot{coda} and scale reduction factors with the R function gelman.diag{coda} using the R package coda [
      • Plummer M.
      • Best N.
      • Cowles K.
      • Vines K.
      CODA: Convergence Diagnosis and Output Analysis for MCMC.
      ]. We are presenting the results as the mean (point estimate) and the 2.5% and 97.5% quantiles (95% credible intervals (CI), [
      • Gelman A.
      • Greenland S.
      Are confidence intervals better termed ǣuncertainty intervalsǥ?.
      ]) of the posterior distribution following [
      • Korner-Nievergelt F.
      • Roth T.
      • Von Felten S.
      • Guélat J.
      • Almasi B.
      • Korner-Nievergelt P.
      Bayesian Data Analysis in Ecology Using Linear Models with R, BUGS, and STAN.
      ].

      Model data set

      The used data set is a European wide gradient study of the impact of climate change and air pollution (N deposition) on forest tree health assessed by foliar nutrient status N:P of different different conifer species (n=10) from the years 1995-2017 [
      • Du E.
      • van Doorn M.
      • de Vries W.
      Spatially divergent trends of nitrogen versus phosphorus limitation across European forests.
      ]. N deposition are based on the EMEP MSC-W model with a 0.1° spatial resolution [
      • Fagerli H.
      • Tsyro S.
      • Jonson J.E.
      • Nyíri Á.
      • Simpson D.
      • Wind P.
      • Benedictow A.
      • Klein H.
      • Mu Q.
      • Denby B.R.
      • Wærsted E.G.
      EMEP - Status Report 1/2020 - Transboundary particulate matter, photo-oxidants, acidifying and eutrophying components.
      ]. Mean annual temperatures (MAT) are based on the E-OBS climate data set with a 0.25° spatial resolution [
      • Du E.
      • van Doorn M.
      • de Vries W.
      Spatially divergent trends of nitrogen versus phosphorus limitation across European forests.
      ].

      Prior settings

      We used the approved critical load for coniferous woodland, which is currently set at 5-15 kg N ha-1 a-1 [
      • Bobbink R.
      • Braun S.
      • Nordin A.
      • Power S.
      • Schütz K.
      • Strengbom J.
      • Weijters M.
      • Tomassen H.
      Review and revision of empirical critical loads and dose-response relationships. Proceedings of an expert workshop, Noordwijkerhout, 23-25 June 2010.
      ], to construct an informative prior (CLempN) assuming a normal distribution with the approved critical load as its mean and half the range as its standard deviation (Normal(mean = 10, sd = 5)) and vague priors with mean = 0 and sd = 2 for fixed effects such as the mean annual temperature (MAT) in this example.

      Bayesian change-point regression model using JAGS (BCR_JAGS)

      According to the current Bayesian framework of change-point regression models by [
      • Roth T.
      • Kohli L.
      • Rihm B.
      • Meier R.
      • Achermann B.
      Using change-point models to estimate empirical critical loads for nitrogen in mountain ecosystems.
      ], we used a linear mixed effect model (c.f. data distribution in Fig. 1) with identity-link function. We described the variation of measured foliar N:P ratio (NPi) between i=1,...,N sampling sites using a Normal distribution with expected foliar N:P ratio λi:
      Fig. 1
      Fig. 1Data distribution of the data set of European foliar N:P ratio. A: Frequency plot of all data including 10 conifer tree species. B: Density distribution grouped by conifer species. C: Density distribution grouped by the different genus types of conifer species, which was used in the final change-point regression model as random effect.
      NPiNormal(λi) and the expected foliar N:P ratio λi as identity link function expressed as:
      λi=β0+k=1KβkXi,k+αi with β0 as intercept and βk as linear slopes for k=1,...,K confounding variables at sampling site i with covariate value Xi,k [
      • Gelman A.
      • Hill J.
      Data analysis using regression and multilevel/hierarchical models.
      ]. According to [
      • Roth T.
      • Kohli L.
      • Rihm B.
      • Meier R.
      • Achermann B.
      Using change-point models to estimate empirical critical loads for nitrogen in mountain ecosystems.
      ] we added αi as the effect of N deposition Ni on the expected foliar N:P ratio λi:
      αi={ 0ifNi<CLβN(NiCL)ifNiCL assuming no effect of N deposition if N deposition Ni is below the critical load (Cl) and a linear change (βN) of λi with increasing N deposition if Ni is equal or above the critical load (CL).
      The random effect term (lines 20-22 in the JAGS model: JAGS_change_point_model_random_effect.R) was defined as:
      λi=β0,s+k=1KβkXi,k+αi
      β0,snorm(μsd) with the intercept β0,s that varies between conifer tree species in this example. The differences between the tree species is modelled using a normal distribution, which is the definition of a random effect.

      BCR_JAGS model diagnostics

      The convergence diagnostics using Gelman and Rubin’s [
      • Brooks S.P.
      • Gelman A.
      General Methods for Monitoring Convergence of Iterative Simulations.
      ] trace and marginal density plots (Fig. 2) showed a good convergence of the two chains. Also the smoothed density plots showed rather balanced histograms of the trace plot values for this model.
      Fig. 2
      Fig. 2Bayesian change-point regression model convergence diagnostics of the model BCR_JAGS, including the proposed integration of the random effect. Left: Trace plots showing the values of the three model parameters during the iterations of the two chains (black and red). Right: Density plot of the three model parameters, showing the distribution of the values in the chains.

      Estimated change-point using JAGS

      The effect plot of the estimated change-point including the random effect of conifer tree species is shown in Fig. 3. The change-point, using the underlying data on foliar N:P ratio and the Bayesian change-point regression model BCR_JAGS, is 5.7±4.8 (Tab.1).
      Fig. 3
      Fig. 3Change-point regression (model BCR_JAGS) of foliar N:P ratios including the fixed effect MAT and tree species as random effect. The grey line represents the estimated critical load () of N including the corresponding 95% credible intervals as dotted lines. Points are measurements from
      [
      • Du E.
      • van Doorn M.
      • de Vries W.
      Spatially divergent trends of nitrogen versus phosphorus limitation across European forests.
      ]
      and are coloured according to conifer tree species. The black line is the estimated change-point regression including 95% credible intervals.

      Bayesian change-point regression model using brms (BCR_brms)

      With the package “brms” [
      • Bürkner P.-C.
      brms : An R Package for Bayesian Multilevel Models Using Stan.
      ] we fitted a Bayesian non-linear multivariate multilevel model with random intercept (β0) and different CLempN between the tree species:
      brmsformula(dependvarβ0+β1*fixeff+step(NdepCLempN)*betaN*(NdepCLempN),CLempN1+(1|random),beta01+(1|random),beta1+betaN1,nl=TRUE,family=gaussian)
      With the same priors and model settings as before in the model BCR_JAGS (chains = 2, iter=100000, thin = 2, warmup = 50000, cores = 4).

      BCR_brms model diagnostics

      The MCMC diagnostics showed a good convergence of the two chains. The posterior distributions are centred around one peak value (Fig. 4).
      Fig. 4
      Fig. 4Bayesian change-point regression model convergence diagnostics of the model BCR_brms. Right: Trace plots showing the two chains. Left: Density plot of the posterior distributions.

      Estimated change-point using brms

      The estimated change-point using the model BCR_JAGS is 9.0±4.3 (Table 1) with an estimated Bayesian R2 = 0.42 [
      • Gelman A.
      • Goodrich B.
      • Gabry J.
      • Vehtari A.
      R-squared for Bayesian Regression Models.
      ]. The effect plot of the change-point is shown in Fig. 5.
      Table 1Estimated critical empirical loads of N with the change-point regression model BCR_JAGS and BCR_brms. Estimated model parameters are given as the mean of the Bayesian posterior distribution and the 95% credible intervals (CI).
      Estimated CLempN (kg ha-1a-1)MeanSDCI 2.5%CI 97.5%
      BCR_JAGS5.714.81-4.6914.81
      BCR_brms9.064.310.4917.76
      Fig. 5
      Fig. 5Change-point regression (model BCR_brms) of foliar N:P ratios including the confounding factor mean annual temperature (MAT) and tree species as random effect. The grey line represents the estimated critical load () of N including the corresponding 95% credible intervals as dotted lines. Points are measurements from
      [
      • Du E.
      • van Doorn M.
      • de Vries W.
      Spatially divergent trends of nitrogen versus phosphorus limitation across European forests.
      ]
      and are coloured according to conifer tree species. The black line is the estimated change-point regression including 95% credible intervals.

      Method validation

      The two presented ways of modelling change-point regressions including random effects, with JAGS (BCR_JAGS) and with the brm function from the brms [
      • Bürkner P.-C.
      brms : An R Package for Bayesian Multilevel Models Using Stan.
      ] package (BCR_brms), resulted in similar estimated CLempN (Table 1), latter being slightly higher with the model (BCR_brms). The estimated CL from the model (BCR_brms) including the non-linearity and the random intercepts of β0 made the outcome more realistic compared to the underlying data (Fig. 5).
      In addition, a simple boxplot comparison (pinpoint method applied by [
      • Bobbink R.
      • Braun S.
      • Nordin A.
      • Power S.
      • Schütz K.
      • Strengbom J.
      • Weijters M.
      • Tomassen H.
      Review and revision of empirical critical loads and dose-response relationships. Proceedings of an expert workshop, Noordwijkerhout, 23-25 June 2010.
      ]) of different groups of N deposition is suggesting a change in foliar N:P ratio between the first group of 0-5 kg ha-1 a-1 and the second group of 5-10 kg ha-1 a-1 (Fig. 6), which verifies the estimated CLempN of both models.
      Fig. 6
      Fig. 6Changes in foliar N:P ratio with increasing N deposition values. The group-wise boxplot comparison (pinpoint method according to
      [
      • Bobbink R.
      • Braun S.
      • Nordin A.
      • Power S.
      • Schütz K.
      • Strengbom J.
      • Weijters M.
      • Tomassen H.
      Review and revision of empirical critical loads and dose-response relationships. Proceedings of an expert workshop, Noordwijkerhout, 23-25 June 2010.
      ]
      ) highlights a suggested change in foliar N:P ratio between the first (0-5 kg ha-1 a-1) and the second group (5-10 kg ha-1 a-1).

      Additional information

      Introduction to critical empirical loads and change-point regression models

      The deposition of atmospheric nitrogen (N) is a major threat to biodiversity and ecosystem functioning globally [
      • Sala O.E.
      • Stuart Chapin F.
      • III
      • Armesto J.J.
      • Berlow E.
      • Bloomfield J.
      • Dirzo R.
      • Huber-Sanwald E.
      • Huenneke L.F.
      • Jackson R.B.
      • Kinzig A.
      • Leemans R.
      • Lodge D.M.
      • Mooney H.A.
      • Oesterheld M.
      • Poff N.L.
      • Sykes M.T.
      • Walker B.H.
      • Walker M.
      • Wall D.H.
      Global Biodiversity Scenarios for the Year 2100.
      ,
      • Stevens C.J.
      How long do ecosystems take to recover from atmospheric nitrogen deposition?.
      ,
      • Bobbink R.
      • Hicks K.
      • Galloway J.
      • Spranger T.
      • Alkemade R.
      • Ashmore M.
      • Bustamante M.
      • Cinderby S.
      • Davidson E.
      • Dentener F.
      • Emmett B.
      • Erisman J.-W.
      • Fenn M.
      • Gilliam F.
      • Nordin A.
      • Pardo L.
      • De Vries W.
      Global assessment of nitrogen deposition effects on terrestrial plant diversity: a synthesis.
      ,
      • Stevens C.J.
      • David T.I.
      • Storkey J.
      Atmospheric nitrogen deposition in terrestrial ecosystems: Its impact on plant communities and consequences across trophic levels.
      ,
      • IPBES
      Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services.
      ]. The accumulation of N over time is one of the main drivers of changes in species distribution across many habitats, especially for habitats with low management regimes such as natural forests [
      • Bobbink R.
      • Braun S.
      • Nordin A.
      • Power S.
      • Schütz K.
      • Strengbom J.
      • Weijters M.
      • Tomassen H.
      Review and revision of empirical critical loads and dose-response relationships. Proceedings of an expert workshop, Noordwijkerhout, 23-25 June 2010.
      ]. The impacts of N deposition are known for certain terrestrial habitats, but many still remains uncertain [
      • Bobbink R.
      • Braun S.
      • Nordin A.
      • Power S.
      • Schütz K.
      • Strengbom J.
      • Weijters M.
      • Tomassen H.
      Review and revision of empirical critical loads and dose-response relationships. Proceedings of an expert workshop, Noordwijkerhout, 23-25 June 2010.
      ]. Critical loads of N has been developed in the framework of the Convention on Long-range Transboundary Air Pollution [
      • Bobbink R.
      • Braun S.
      • Nordin A.
      • Power S.
      • Schütz K.
      • Strengbom J.
      • Weijters M.
      • Tomassen H.
      Review and revision of empirical critical loads and dose-response relationships. Proceedings of an expert workshop, Noordwijkerhout, 23-25 June 2010.
      ,
      • de Vries W.
      • Posch M.
      • Sverdrup H.U.
      • Larssen T.
      • de Wit H.A.
      • Bobbink R.
      • Hettelingh J.-P.
      Geochemical Indicators for Use in the Computation of Critical Loads and Dynamic Risk Assessments.
      ]. These critical loads are defined as thresholds below which damaging effects on specific habitats do not occur based on the latest scientific knowledge [

      J. Nilsson, P. Grennfelt, N.C. of Ministers, Critical loads for sulphur and nitrogen. Report from a workshop held at Skokloster, Sweden, 19-24 March, 1988, 1988.

      ]. Critical loads of N are often defined based on negative effects on plant diversity [
      • Bobbink R.
      • Hicks K.
      • Galloway J.
      • Spranger T.
      • Alkemade R.
      • Ashmore M.
      • Bustamante M.
      • Cinderby S.
      • Davidson E.
      • Dentener F.
      • Emmett B.
      • Erisman J.-W.
      • Fenn M.
      • Gilliam F.
      • Nordin A.
      • Pardo L.
      • De Vries W.
      Global assessment of nitrogen deposition effects on terrestrial plant diversity: a synthesis.
      ], as they are often lower compared to critical loads of N from soil processes such as soil acidification or N leaching [
      • Du E.
      Effects of Nitrogen Deposition on Forest Ecosystems.
      ]. Empirical critical loads can be determined with the use of N addition or reduction experiments or in gradient studies covering a gradient of N deposition. Latter is a more holistic approach taking into account the variability of natural habitats. However, gradient studies can also have various limitations for instance they should use a sufficient spatial scale of the air pollutant assessed and they should take into account the most important modifying factors [
      • Braun S.
      • Schindler C.
      • Rihm B.
      Growth trends of beech and Norway spruce in Switzerland: The role of nitrogen deposition, ozone, mineral nutrition and climate.
      ]. In addition, an appropriate framework of regression model should be selected, suitable for the selected study design. Bayesian change-point regression models have been shown to overcome some difficulties of large spatial variation found in natural habitats by taking into account confounding factors such as climatic variation or different soil types [
      • de Vries W.
      • Posch M.
      • Sverdrup H.U.
      • Larssen T.
      • de Wit H.A.
      • Bobbink R.
      • Hettelingh J.-P.
      Geochemical Indicators for Use in the Computation of Critical Loads and Dynamic Risk Assessments.
      ,
      • Groffman P.M.
      • Baron J.S.
      • Blett T.
      • Gold A.J.
      • Goodman I.
      • Gunderson L.H.
      • Levinson B.M.
      • Palmer M.A.
      • Paerl H.W.
      • Peterson G.D.
      • Poff N.L.
      • Rejeski D.W.
      • Reynolds J.F.
      • Turner M.G.
      • Weathers K.C.
      • Wiens J.
      Ecological Thresholds: The Key to Successful Environmental Management or an Important Concept with No Practical Application?.
      ], allowing a more accurate estimation of the critical load [
      • Beckage B.
      • Joseph L.
      • Belisle P.
      • Wolfson D.B.
      • Platt W.J.
      Bayesian change-point analyses in ecology.
      ]. Therefore, change-point regression models applied in a Bayesian framework are useful statistical tools in estimating critical empirical loads [
      • Roth T.
      • Kohli L.
      • Rihm B.
      • Meier R.
      • Achermann B.
      Using change-point models to estimate empirical critical loads for nitrogen in mountain ecosystems.
      ].

      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.

      Acknowledgments

      This paper is a collaborative piece of work between ST and SB from the Institute for Applied Plant Biology providing the idea for the adaption of the current change-point method in the frame work of the revision of the critical empirical load for nitrogen and TR from Hintermann & Weber AG providing the statistical back ground. ED provided the data and edited the manuscript. This project was supported by the Swiss Federal Office for the Environment (FOEN).

      Appendix A. Supplementary materials

      References

        • Roth T.
        • Kohli L.
        • Rihm B.
        • Meier R.
        • Achermann B.
        Using change-point models to estimate empirical critical loads for nitrogen in mountain ecosystems.
        Environ. Pollut. 2017; 220: 1480-1487https://doi.org/10.1016/j.envpol.2016.10.083
        • Du E.
        • van Doorn M.
        • de Vries W.
        Spatially divergent trends of nitrogen versus phosphorus limitation across European forests.
        Sci. Total Environ. 2021; 771: 145391https://doi.org/10.1016/j.scitotenv.2021.145391
      1. R Core Team, Team R Development Core, R: A Language and Environment for Statistical Computing, 2021, http://www.r-project.org/.

      2. J.J. Allaire, Y. Xie, J. McPherson, J. Luraschi, K. Ushey, A. Atkins, H. Wickham, J. Cheng, W. Chang, R. Iannone, rmarkdown: Dynamic Documents for R, 2021. https://github.com/rstudio/rmarkdown.

        • Plummer M.
        JAGS: A Program for Analysis of Bayesian Graphical Models using Gibbs Sampling.
        3rd Int. Work. Distrib. Stat. Comput. (DSC 2003); Vienna, Austria. 2003; 124
      3. M. Plummer, rjags: Bayesian Graphical Models using MCMC, 2019. https://cran.r-project.org/package=rjags.

        • Bürkner P.-C.
        brms : An R Package for Bayesian Multilevel Models Using Stan.
        J. Stat. Softw. 2017; 80https://doi.org/10.18637/jss.v080.i01
        • Roth T.
        • Sprau P.
        • Naguib M.
        • Amrhein V.
        Sexually selected signaling in birds: A case for Bayesian change-point analysis of behavioral routines.
        Auk. 2012; 129: 660-669https://doi.org/10.1525/auk.2012.12041
        • Link W.A.
        • Cam E.
        • Nichols J.D.
        • Cooch E.G.
        Of Bugs and Birds: Markov Chain Monte Carlo for Hierarchical Modeling in Wildlife Research.
        J. Wildl. Manage. 2002; 66: 277https://doi.org/10.2307/3803160
        • Brooks S.P.
        • Gelman A.
        General Methods for Monitoring Convergence of Iterative Simulations.
        J. Comput. Graph. Stat. 1998; 7: 434-455https://doi.org/10.1080/10618600.1998.10474787
        • Plummer M.
        • Best N.
        • Cowles K.
        • Vines K.
        CODA: Convergence Diagnosis and Output Analysis for MCMC.
        R News. 2006; 6: 7-11
        • Gelman A.
        • Greenland S.
        Are confidence intervals better termed ǣuncertainty intervalsǥ?.
        BMJ. 2019; 366: l5381https://doi.org/10.1136/bmj.l5381
        • Korner-Nievergelt F.
        • Roth T.
        • Von Felten S.
        • Guélat J.
        • Almasi B.
        • Korner-Nievergelt P.
        Bayesian Data Analysis in Ecology Using Linear Models with R, BUGS, and STAN.
        Elsevier, 2015https://doi.org/10.1016/C2013-0-23227-X
        • Fagerli H.
        • Tsyro S.
        • Jonson J.E.
        • Nyíri Á.
        • Simpson D.
        • Wind P.
        • Benedictow A.
        • Klein H.
        • Mu Q.
        • Denby B.R.
        • Wærsted E.G.
        EMEP - Status Report 1/2020 - Transboundary particulate matter, photo-oxidants, acidifying and eutrophying components.
        Technical Report. Norwegian Meteorological Institute, 2020
        • Gelman A.
        • Hill J.
        Data analysis using regression and multilevel/hierarchical models.
        Cambridge university press, 2006
        • Gelman A.
        • Goodrich B.
        • Gabry J.
        • Vehtari A.
        R-squared for Bayesian Regression Models.
        Am. Stat. 2019; 73: 307-309https://doi.org/10.1080/00031305.2018.1549100
        • Sala O.E.
        • Stuart Chapin F.
        • III
        • Armesto J.J.
        • Berlow E.
        • Bloomfield J.
        • Dirzo R.
        • Huber-Sanwald E.
        • Huenneke L.F.
        • Jackson R.B.
        • Kinzig A.
        • Leemans R.
        • Lodge D.M.
        • Mooney H.A.
        • Oesterheld M.
        • Poff N.L.
        • Sykes M.T.
        • Walker B.H.
        • Walker M.
        • Wall D.H.
        Global Biodiversity Scenarios for the Year 2100.
        Science (80-.). 2000; 287: 1770-1774https://doi.org/10.1126/science.287.5459.1770
        • Stevens C.J.
        How long do ecosystems take to recover from atmospheric nitrogen deposition?.
        Biol. Conserv. 2016; 200: 160-167https://doi.org/10.1016/j.biocon.2016.06.005
        • Bobbink R.
        • Hicks K.
        • Galloway J.
        • Spranger T.
        • Alkemade R.
        • Ashmore M.
        • Bustamante M.
        • Cinderby S.
        • Davidson E.
        • Dentener F.
        • Emmett B.
        • Erisman J.-W.
        • Fenn M.
        • Gilliam F.
        • Nordin A.
        • Pardo L.
        • De Vries W.
        Global assessment of nitrogen deposition effects on terrestrial plant diversity: a synthesis.
        Ecol. Appl. 2010; 20: 30-59https://doi.org/10.1890/08-1140.1
        • Stevens C.J.
        • David T.I.
        • Storkey J.
        Atmospheric nitrogen deposition in terrestrial ecosystems: Its impact on plant communities and consequences across trophic levels.
        Funct. Ecol. 2018; 32: 1757-1769https://doi.org/10.1111/1365-2435.13063
      4. Bonn, Germany
        • de Vries W.
        • Posch M.
        • Sverdrup H.U.
        • Larssen T.
        • de Wit H.A.
        • Bobbink R.
        • Hettelingh J.-P.
        Geochemical Indicators for Use in the Computation of Critical Loads and Dynamic Risk Assessments.
        Crit. Loads Dyn. Risk Assessments. 2015: 15-58https://doi.org/10.1007/978-94-017-9508-1_2
      5. J. Nilsson, P. Grennfelt, N.C. of Ministers, Critical loads for sulphur and nitrogen. Report from a workshop held at Skokloster, Sweden, 19-24 March, 1988, 1988.

        • Du E.
        Effects of Nitrogen Deposition on Forest Ecosystems.
        Handb. Air Qual. Clim. Chang. Springer Singapore, Singapore2022: 1-23https://doi.org/10.1007/978-981-15-2527-8_27-1
        • Braun S.
        • Schindler C.
        • Rihm B.
        Growth trends of beech and Norway spruce in Switzerland: The role of nitrogen deposition, ozone, mineral nutrition and climate.
        Sci. Total Environ. 2017; 599-600: 637-646https://doi.org/10.1016/j.scitotenv.2017.04.230
        • Groffman P.M.
        • Baron J.S.
        • Blett T.
        • Gold A.J.
        • Goodman I.
        • Gunderson L.H.
        • Levinson B.M.
        • Palmer M.A.
        • Paerl H.W.
        • Peterson G.D.
        • Poff N.L.
        • Rejeski D.W.
        • Reynolds J.F.
        • Turner M.G.
        • Weathers K.C.
        • Wiens J.
        Ecological Thresholds: The Key to Successful Environmental Management or an Important Concept with No Practical Application?.
        Ecosystems. 2006; 9: 1-13https://doi.org/10.1007/s10021-003-0142-z
        • Beckage B.
        • Joseph L.
        • Belisle P.
        • Wolfson D.B.
        • Platt W.J.
        Bayesian change-point analyses in ecology.
        New Phytol. 2007; 174: 456-467https://doi.org/10.1111/j.1469-8137.2007.01991.x