Covering an area of 8,700 km², Crete is the largest and most populated island of Greece, and the fifth largest in the Mediterranean. According the Köppen classification, Crete has a Mediterranean—Semiarid climate featuring long and dry summers, and relatively wet and cold winters (Kottek et al., 2006). Crete receives on average about 7.7 billion m³ of rainfall, of which only and 10–15% produces runoff, while 68–76% evapotranspires and 14–17% infiltrates (Koutroulis et al., 2016). The intense tectonic history has formed the island's complex topography that ranges from sea level to 2,450 m, and is abundant in small, ephemeral watersheds (Tsanis et al., 2011). This highly-rugged terrain has been definitive for human development as well as its spatial allocation (Koutroulis et al., 2016). Similarly, this variability has spurred the development of the wide variety of biotopes present on the island (Sfikas, 1987), ranging from coastal to alpine, and the reciprocal plethora of endemic and rare species that constitute one of the 10 biodiversity hotspots in the Mediterranean (Medail and Quezel, 1997) numbering 1,624 native and 47 introduced vascular floral species (Turland et al., 1993). The importance of this biodiversity is highlighted by the fact that over 30% of the island has been included in the Natura 2000 protected area network (Dimitrakopoulos et al., 2004).

Maps of presence-absence of vascular plant species distributions in Crete were digitized from Turland et al. (1993) and its latest supplement (Chilton and Turland, 2004). The island of Crete and its surrounding islets were divided into 162 grid cells, each covering an area of 8.25 × 8.25 km², following the grid cell size of Turland et al. (1993). On each cell, the native, endemic, and alien species richness was calculated. We used (Turland et al., 1993; Chilton and Turland, 2004; and references therein) to define native (nnat = 1,395) and endemic (nend = 174) species, and the vascular plants from D'Agata et al. (2009) that are listed in Chilton and Turland (2004) and Turland et al. (1993) were used to define alien species richness. Only species present in at least two cells were used (nalien = 47). Coarse-grid information was estimated by aggregating this dataset by a factor of two, thus reducing resolution to grid cells of 16.5 × 16.5 km². The spatial distribution of the original as well as the resampled data regarding alien species presences are visualized in Figure 1. All input variables and their ranges (min – max values within each cell) are listed in Table 1.

Habitat classification relied on the most detailed resolution available of the CORINE Landcover (level 3, spatial resolution 100 m; EEA-ETC/TE., 2002), to calculate the richness and percentage of every land cover class within every grid cell, using Patch Analyst 5.1 within ArcGIS. In order to avoid potential temporal deviance between habitat classifications and species presences in cells, the last updated available supplement for the flora of Crete published in 2008 (Chilton and Turland, 2004) and the closest available time snapshot of the CORINE landcover for Crete in 2010 were used. The classification process resulted in 29 habitat types, of which 9 agricultural, 7 artificial, and 13 natural. We recorded habitat richness per cell as the number of different land cover types present on each cell (total, artificial, agricultural, and natural habitat richness) as well as percentage of cover (total, artificial, agricultural, and natural % of cell cover).

Climatic variables were derived from WorldClim (Hijmans et al., 2005) for Crete and surrounding islets. The original resolution of the climatic data was 1 × 1 km². In order to re-scale them to 8.25 km and match them with the grid of the plant data, the mean values of the 1 km data within the 8.25 km cells were calculated and used. The climatic variables used were annual mean temperature (Tempmean), annual mean temperature of warmest quarter (Tempwarm), annual mean temperature of coldest quarter (Tempcold), all in °C, annual mean precipitation (Precipmean), precipitation of wettest quarter (Precipwet), and precipitation of driest quarter (Precipdry), all in mm year⁻¹. Soil data were derived from SoilGrid (Hengl et al., 2014) and rescaled from 1 to 8.25 km as the climatic data. The soil variable used was soil richness in the cell (Soildiv) derived as the number of different soil types occurring within each cell. The indices of elevation recorded were the mean of all elevation values within the cell (Alt) and the range of elevation within the cell (Alt range) both in meters.

Random Forests (RFs; Breiman, 2001) take advantage of boosting (Schapire et al., 1998) and bagging (bootstrap aggregating; Breiman, 1996a) of the Classification And Regression Tree (CART; Breiman et al., 1984) model, and adapt a more random but nevertheless more efficient node splitting strategy than standard CARTs (Liaw and Wiener, 2002). In RFs, each individual tree is developed after the following steps: (1) Given a set of training data N, n random samples with repetition (bootstrap) are taken as training set; (2) For each node of the tree, M input variables are determined, and m ≪ M, variables are selected for each node. The most important variable randomly chosen is used as a node. The value of m remains constant; (3) Each tree is developed to its maximum expansion.

RFs have been employed in a wide variety of classification and prediction problems (Scornet et al., 2015; Cano et al., 2017) as they are among the most effective computationally-intensive algorithms to extract information from unstable estimates (Scornet et al., 2015). They are especially well suited for large, high-dimensional datasets, where problem complexity and scale render direct discovery of a good model in a single step impossible (Büchlmann and Yu, 2002; Kleiner et al., 2014; Wager et al., 2014). The fact that RFs require tuning of only two parameters (the tree population in each forest and the number of input variables m randomly selected at each node) for which they are usually not very sensitive (Liaw and Wiener, 2002), and their accuracy and competence when faced with scarce, multivariate datasets of intricate structure (Scornet et al., 2015), have greatly contributed to their popularity.

Similar to other data-driven approaches, RFs may not perform equally well when the task at hand is extrapolatory beyond the range of the recovered predictor-predictand relationship or involves scenario analysis (Daliakopoulos and Tsanis, 2016). Furthermore, Strobl and Boulesteix (2007) showed that variable importance measures of the original RF algorithm may be biased due to differences among predictor structure and scale, adding to the interpretability challenges of data-driven methods. Nevertheless, an extensive data-driven model inter-comparison by Fernández-Delgado et al. (2014) showed that they may be the first weapon of choice for real-world problems.

Typically, CARTs error is estimated following the out-of-bag (OOB) error R(D) of a selection of the input observations based on bagging, otherwise an OOB sample D (James et al., 2013). In RFs, for each tree t, prediction error of D is estimated before and after randomly permuting the values of the j-th variable, thus giving R(Dnt) and R(Dntj), respectively. Typically, imbalanced datasets favor correct classification of the majority class, nevertheless, RFs can account for this bias by adjusting the voting cut-off from the default 1/c, where c is the number of classes. This provides additional flexibility to the RF algorithm (Ma et al., 2006) and allows for favoring sensitivity or specificity to different classes. A variable can be considered a strong predictor when permuting it increases the prediction error (Gregorutti et al., 2017), therefore it's importance I_V can be defined as:

The Mean Decrease in Accuracy (MDA) is estimated by averaging this difference over all trees, and normalizing it by the standard deviation of the differences. The more the accuracy of the RF decreases due to the exclusion (or permutation) of a single predictor, the more important that predictor is considered, and therefore variables with a large MDA are more important for data classification.

Gini is one of the most encountered impurity functions, providing a measure of the “goodness-of-split” for CARTs by favoring splits that allocate a single pure node for the largest class and the rest for the remaining classes (Breiman, 1996b). The Gini index for a node t can be calculated as:

where c is the number of classes and p(i|t), p(j|t) are the estimated probabilities of classes i, j at node t (Cano et al., 2017). In this context, Mean Decrease Gini (MDG) aggregates the Gini gain over all splits and trees to assess the classifying capacity of a variable (Friedman et al., 2009) and is thus a metric of the homogeneity of nodes and leaves in the RF (Bluemke and Stepień, 2016).

MDA and MDG can rank each independent variable for its effectiveness as a predictor of alien species richness, but don't show or quantify the actual positive, negative, humped, etc. relationship between them. Nevertheless, this is an elementary process under conditions of multiple acting variables (Häring et al., 2012), such as cumulative human impacts. For this reason, partial dependence plots (Friedman, 2001; Friedman et al., 2009) can be used to depict the relationship of alien species presence probability on each predictor after averaging out the effects of all classification predictors (Cutler et al., 2007).

Finally, the Receiver Operating Characteristics (ROC) analysis has been an indispensable tool for signal detection and diagnostic systems. As documented by Pontius and Si (2014), ROC has been employed in a wide range of applications requiring a threshold-independent measure to compare predicted against observed values. ROC plots have been previously considered in plant ecology, both at a theoretical (Guisan and Zimmermann, 2000) and applied (Manel et al., 2002; Wang et al., 2014) level as effective indicators of model performance independent of the threshold probability. Typically, ROC curves depict true positive rate (TPR), otherwise sensitivity, recall or hit rate, against true negative rate (TNR), otherwise called specificity. In terms of model estimates, TPR and TNR are defined as:

where T, F, P, and N stand for true, false, positive and negative, respectively. The complementary values of TPR and TNR are false negative rate (FNR), otherwise miss rate, and false positive rate (FPR), otherwise fall-out or false alarm. Based on these values, the Matthew's correlation coefficient (MCC; Matthews, 1975), a reduction of the Pearson correlation coefficient for binary variables (Baldi and Brunak, 2001), is a popular evaluation criterion of machine learning performance (Bhasin and Raghava, 2004; Chen et al., 2004; Bao and Cui, 2005):

MCC has an advantage in imbalanced datasets where the disparity in the number of presence and absence samples is significant.

Experiments were developed using the latest (v4.6–12) implementation of Breiman and Cutler's original Fortran code by Liaw and Wiener (2002) in R. While RFs can be trained very efficiently and avoid overfitting (Breiman, 2001), predictions and variable significance ranking are seldom the identical after each random training, especially for small datasets. To account for this uncertainty, a bootstrapping approach of training multiple RFs is adopted. For each training iteration k, RF_k is presented with 70% of the dataset, sampled with replacement, and the remaining is reserved for testing. Presenting only part of the dataset to the RFs also simulates operational use where only part of the study area is sampled at fine grid and the rest is sampled at coarse-grid resolution. Furthermore, as subsets of alien species presence and absence were imbalanced, training was executed using a variable training cutoff, ranging from 0.1 to 0.9. The full code in R used for the analysis is provided in Supplementary Material.

Mean decrease in accuracy (MDA) results as estimated from bootstrap randomizations indicate that, apart from the coarse resolution alien species presence, the percentage of natural cover within each cell was the most important predictor of alien species presence, followed by the endemic species richness, altitude, minimum temperature, and altitude range within each cell (Figure 2, left). From the ones explored here, the least predictive in MDA were artificial habitat richness, temperature range, habitat richness, the percentage of the surface area of each cell within the Natura 2,000 protected area network, and the soil type richness (Figure 2, left). In the latter cases, some bootstrap samples have yielded negative results suggesting that permuting these variables from the predictor vector increases accuracy. Results in Mean Decrease Gini (MDG) are in general agreement with those of MDA, also evaluating natural cover, endemic species richness, and altitude as the most efficient splitting variables (Figure 2, right). Agricultural cover replaces minimum temperature for the MDG rating but both variables score highly for both criteria. The least efficient node splits according to MDG were performed by artificial habitat richness, natural habitat richness, agricultural habitat richness, soil type richness, and temperature range (Figure 2, right). Emphatically, artificial habitat richness is the worst predictor for both metrics, essentially boosting the noise in the dataset.

Results from partial dependence plots among the most predictive variables according to the MDA and MDG criteria indicate that the percentage of natural cover has an overall positive relationship with alien species richness, while the percentage of agricultural cover has an overall negative relationship with alien species richness (Figure 3). Altitude, and altitude range has an overall positive relationship between alien species richness, mean temperature has a negative relationship for larger temperature values while mean annual precipitation and precipitation range has a humped relationship with alien species richness (Figure 3). Therefore, for the case at hand, in the event of a survey priority may be given to low-temperature, elevated natural areas with high topographic variability, far from agricultural use and precipitation extremes. The percentage of each cell within the Natura 2,000 protected area network has a positive relationship with alien species richness (Figure 3), albeit this variable was not within the most predictive of alien species richness based on MDA or MDG.

As shown by the results, the bootstrapping method followed herein is helpful for drawing a more robust conclusion, particularly regarding the partial dependence plots. Bootstrapped predictors (solid black lines in Figure 3) are more stable, less prone to overfit and more inclusive than single experiment predictors. This becomes obvious in the Mean Temperature plot of Figure 3, where the red line representing an OOB sample does not cover the entire range of temperature values in the dataset. As low temperatures are not common in the dataset, the OOB estimation of dependence does not always include these values. Using an additional layer of bootstrapping ensures that the full range of values is explored.

True negative detection rates (TNR; not detecting alien species in cells where alien species are not present) declines with an increasing cut-off rate while true positive detection rates (TPR; detecting alien species in cells where alien species are present) increases with an increasing cut-off rate (Figure 4). MCC values indicate a strong positive relationship at cut-offs between 0.2 and 0.5 and are otherwise acceptable correlation. When cut-off increases TNs decrease and TFs increase, therefore more alien species can be detected but by being more exhaustive more false alarms are also generated. When cut-off remains low, less risk is taken with surveying resources but a significant fraction of alien species presences is missed.

The author ID was employed by company TM Solutions. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.