You used a two-step process for variable selection, first using AIC for MODIS buffer selection and then VIF for bioclimatic variables after a preliminary Random Forest run. This sequential, species-specific selection approach carries a high risk of overfitting and data dredging, as it maximizes model fit for the training data. More critically, this process may remove variables that are important for understanding nonstationarity (e.g., bioclim_18 was important for NOFL in NWFM but not MWCF). How do you justify this selective process given your core aim is to compare relative variable importance across ecoregions? The differing variable sets complicate direct comparisons.
For rare species (BBWO, WHWO), your effective sample sizes are very small and detection frequencies are ~1%. The training dataset sizes for these species are also small, and they are heavily influenced by the spatial subsampling (e.g., BBWO went from 915 to 560 detection records). With such sparse data, your Random Forest models are highly unstable. The high specificity (1.00) and low sensitivity (0.21) for BBWO suggest the model is simply predicting absence. How confident are you that these models are capturing genuine ecological relationships rather than noise or spatial artifacts from the subsampling?
You consistently found that bioclimatic variables (coarse resolution) outranked GEDI-fusion variables (finer resolution) at home-range scales. While you suggest this is due to lower-order habitat selection, it could be a statistical artifact: the coarse bioclimatic variables (800 m) capture broad environmental gradients that are less variable within a species’ home range compared to the finer-scaled GEDI metrics, making them inherently stronger predictors at this scale. Did you test for multicollinearity between the selected bioclimatic variables and the GEDI metrics? If not, how can you be sure the higher importance of bioclimatic variables isn’t simply due to them acting as proxies for the structural variables?