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).