Environmental predictors
To estimate the limiting effect of environmental gradients on species maximum abundances, we used growing degree days (GDD – °C) to measure energy available; and the difference between precipitation and potential evapotranspiration (P-PET; Zang et al., 2020) as a proxy for water balance. To calculate these variables, we obtained temperature and precipitation data in raster format from the Global Meteorological Forcing Dataset for land surface modeling (Sheffield et al., 2006, available in http://hydrology.princeton.edu/data/pgf/v3) spanning from 1965-2016. For details on how we processed the environmental predictors, see Appendix S3.
Quantile regressions
We ran QR specified at upper quartiles to estimate the influence of GDD and water balance on the maximum abundance of each species. Parameterτ was defined for each species to be proportional to its prevalence across the study region (i.e., number of ground plots with positive abundance, so that the greater the number of observations, the better the estimate of the upper limit will be), so that\(\tau=\frac{40}{\text{prevalence}}\) (Villén-Pérez et al. , 2022). Thus, τ parameter varied from 0.75 to 0.998 (Appendix S4 and S5). To account for possible non-linear relationships, we defined both linear and quadratic terms for the environmental predictors. Prior to the calculation of quadratic terms, we z-standardized the original variables, transforming their absolute values to a scale in which the sample mean is 0 and the standard deviation is 1. This z-standardization allows us to compare the relative contribution of each independent predictor in the estimation of the response variable (see Carrascal et al., 2016 for a similar approach), as well as to reduce the problem of multicollinearity between linear and quadratic terms (Kim, 1999). To avoid negative abundance predictions, we summed a small number to all observed abundance values (half of the smallest observed non-0 abundance) and logit-transform them prior to applying QR. Then, we back-transformed the estimated values and subtracted the small number (Orsini & Bottai, 2011; Villén‐Peréz et al. , 2020). The analyses were performed in R (R Core Team, 2019) using package quantreg(Koenker, 2019). After fitting models, we estimated Variation Inflation Factors (VIF; Belsey, 1991) to assess the effect of multicollinearity and only accepted models with VIF < 10 (Neter, 1996) for each of the predictors (Appendix S4 and S5). We also compared each model with the corresponding null model through Akaike information criterion (AIC) and calculated the Evidence Ratio (ER), a measure of how much more likely models including predictors were better than null models (i.e., intercept-only models; Burnham & Anderson, 2002; Symonds & Moussalli, 2011).