The evaluation of the models showed that both AUC and adjusted explained deviance values increased with larger background radius (Fig. 2). 34 models were above the selection thresholds for both AUC (0.65) and adjusted explained deviance (0.15). The average AUC values for the best models was 0.76 (0.05 SD), while average explained deviance was 0.21 (0.05 SD).
The final habitat suitability model (Fig. 3) showed high habitat suitability values and low standard deviation (high agreement among models) in central and southern France, Italian Peninsula, the Mediterranean Islands, along the Atlantic coast of the Iberian Peninsula, as well as along the southern coast of the Black Sea and the Levantine coasts of the Mediterranean Sea. Medium habitat suitability values were found at the continental areas of the Iberian Peninsula, British islands, Dinaric Alps, Balcanic Mountains, continental areas of Anatolia and the southern coasts of the Baltic Sea. Very low habitat suitability values consistent with non-suitable habitats were dominant across the high-altitude regions of the Alps and Pyrenees, across the Scandinavian Shield, and the current Arabian Desert. The lack of agreement among the best models was high in the north of the British Isles, the Scandinavian uplands, the margins of the Alps, and along the northern limits of the Arabian desert.
The habitat suitability values of the presence records (Fig. 4) was generally high (average 0.69, SD 0.14), but most sites at the north and eastern distribution limits (Germany, Slovakia, Hungary, Ukraine, Turkey) showed low habitat suitability values. The average standard deviation for the Neanderthal sites was 0.12 (SD 0.03), but two sites from Slovakia showed higher deviation values: Gánovce (GANOV) and Hôrka-Ondrej (HORK).
The analysis of variable importance performed with Random Forest was robust (95.23 explained deviance), and showed that temperature of the warmest month (69.82 % increment in mean squared error), annual rainfall (57.28 %IncMSE), and temperature of the coldest month (56.69 %IncMSE) were the most important factors shaping habitat suitability at the continental scale. Topographical slope (34.10 %IncMSE), topographic diversity (32.46 %IncMSE) and precipitation of the warmest month (27.17 %IncMSE) were the least important variables. The response curves of the most important predictors (Fig. 5) showed that optimum habitat suitability values occurred when the temperature of the warmest month was between 25 and 39 ºC, the annual rainfall was higher than 700 mm, and the temperature of the coldest month was higher than -5ºC.
The recursive partition analysis of the 44 localities revealed the different environmental processes shaping habitat suitability across Europe. The model based on variable values (Fig. 6) showed that a minimum temperature of the coldest month around -4.4 ºC separated localities with low and high habitat suitability. Localities within the suitable conditions and a topographical slope higher than 2.5º showed the highest habitat suitability, while localities with lower topographical slopes showed intermediate suitability values.
The model based on the local R-squared of the predictors (Fig. 7) showed terminal nodes with a wide range of habitat suitability values, and many cases with a high residual (high difference between actual and predicted suitability according to the recursive partition model). Node 2 grouped low habitat suitability localities in which the temperature of the warmest month was the key factor, correspond to areas either hot or cold during the summer, like the Pyrenees (ID=2) or Central Anatolia (ID=5). The split between localities with moderate and high habitat suitability values was defined by the importance of topographic diversity.
Finally, the model based on the local coefficients of the predictors (Fig. 8) showed that a positive coefficients of temperature of the coldest month (interpreted as sites not warm enough) grouped localities with low habitat suitability (node 2), but a few cases showed low reliability due to their high habitat suitability values (Cantabric Coast, ID=25; Corsica, ID=23). Higher habitat suitability values (node 7) were defined by negative or neutral regression coefficients for the minimum temperature of the coldest month, the topographical slope and the precipitation of the warmest month, indicating in every case variable values around the optimum, or slightly beyond the optimum.