The study area (20ºN to 70ºN, 10ºW to 70ºE) comprises Europe and the Irano-Turanian region. It does not cover the entire known range of Neanderthals, which extended at least to Okladnikov in Southern Siberia \cite{Krause_2007}, but there is only one single record (Ust-Izhul, Siberia) attributed to MIS 5e beyond the Caspian Sea \cite{Rolland2010}, and therefore the data available would be too scarce to build a reliable model for this region.
We reviewed the literature searching for archaeological sites with Neanderthal fossil or lithic remains attributed to MIS 5e, and found 50 records located in fourteen countries, 25 of them with absolute dates, and 25 of them with relative dates (see Table 1). To reduce spatial clustering we filtered the data applying a minimum distance of 100 km between nearby locations \cite{Guisan_2005}, which reduced the sample size to 29 presence records (see Fig. 1).
We obtained the climatic predictors from the Eemian (\(\sim\)130 ka BP) palaeoclimatic simulation produced by \citet{Otto2006}. The data is available at www.worldclim.com \cite{Hijmans20051965} as a set of maps representing 19 bioclimatic variables downscaled to 1km resolution. We derived maps of slope, topographic wetness index, and topographic diversity from the digital elevation model available in the same source. All the variables were aggregated to 5 km resolution, and the sea level was set at 7 m.a.s.l \cite{Dutton_2012}. We assessed multicollinearity among the 19 bioclimatic and 3 topographic predictors using a variance inflation factor analysis (R package HH, \cite{heiberger2015hh}). We iteratively removed correlated variables with the objective of retaining predictors related with extreme temperatures, water availability and topography. The final set of uncorrelated predictors comprised: maximum temperature of warmest month (bio5), minimum temperature of coldest month (bio6), annual precipitation (bio12), precipitation of warmest quarter (bio18), slope and topographic diversity.
Observing all potential biases in our data (lack of absences, low number of presences, uncertainty on the dates of the localities, unknown sampling bias), we selected generalized linear models (GLMs) calibrated with weighted background \cite{Barbet_Massin_2012} and second degree polynomials as modeling method. The use of a weighted background compensates for the lack of absences by providing a comprehensive sampling on the ecological conditions across the study area. However, the geographical distribution of the background points affect modeling outcomes \cite{VanDerWal_2009,Chefaoui_2008}, and a few studies recommend to restrict the background to areas within the dispersal distance of the focus species \cite{Barve_2011,Acevedo_2012}. In our case, there was no available data to develop such criteria, and we therefore explored the sensitivity of the model to different sets of background points generated within buffers around the presence records at increasing radius (100, 200, 400, 800 and 1600 km). Each background dataset sampled the 20% of all the cells within each buffer.
To avoid over-parameterization when fitting a GLM, at least ten presence points per predictor are required \cite{opac-b1127929}. Considering our data (six predictors and 29 presence records), the maximum number of predictors we could use to fit a model was just 3, so there was a high risk of fitting a model missing important predictors. To overcome this problem, we used the dredge function of the R package MuMIn \cite{librarymumin} to generate all the GLM equations combining the six predictors in groups of one, two and three, resulting in 41 different equations. We calibrated one model for each combination of equation and background radius, resulting in 205 different models to evaluate.
We applied a leave-one-out approach to compute AUC values based on 1000 pseudoabsences not used to calibrate the model. This allowed to evaluate relative omission error \cite{Fielding199738,Phillips_2008}, and adjusted explained deviance to assess goodness of fit \cite{Guisan_1999}. We evaluated a total of 5330 models: 41 equations per 5 background radiuses per 26 presence records. We aggregated the AUC and adjusted explained deviance values by equation and background radius. All models with an average AUC higher than 0.65 and average adjusted explained deviance higher than 0.15 were calibrated again using the complete dataset. We computed the average (habitat suitability) and the standard deviation (measure of model agreement) of the habitat suitability maps resulting from these models to produce the final ensemble \cite{ARAUJO_2007,Marmion_2009}. We mapped the average and the standard deviation into a single map using the whitening method proposed by \citet{Hengl200475} to enhance visualization and model interpretation.
Random Forest, a machine learning method based on the ensemble of classification or regression trees \cite{Breiman20015}, is regarded as a robust method to assess variable importance \cite{Cutler20072783}. We applied the randomForest R function \cite{randomforestcitation} to assess the influence of the environmental factors over habitat suitability at the continental scale. This analysis also provided a set of response curves produced by the plotmo R library, \citet{plotmocitation}), useful to understand the effect of each predictor over habitat suitability.
The analyses described in this paper so far cannot evaluate the importance of a particular variable at local scales. This is relevant as habitat suitability might be determined by different predictors at smaller spatial scales. Thus, we define local scale as the average home range of Neanderthals which, according to Daujeard et al. (2012), and based on the transportation of raw lithic materials, was around 50 kilometers \cite{FéblotAugustins1993211}. In consequence we divided the study area into cells of 50 per 50 kilometers, and for each cell with more than 30 original cells (the resolution of the habitat suitability map and the predictors) we fitted one linear regression model per environmental predictor using habitat suitability as response variable. For any given predictor, we considered the R squared to be an indicator of its importance at the local scale, and the coefficient to be an indicator of its effect over habitat suitability. We interpreted near zero coefficients linked to low habitat suitability as regions with extreme values for the given predictor, while near zero coefficients linked to high habitat suitability values were interpreted as optimum habitat. Positive coefficients indicated a positive (but still sub-optimum) effect of the given predictor over habitat suitability, while negative coefficients indicated that the predictor values were beyond the optimum (e.g. too hot, too wet). All the results with p-values higher than 0.05 were recoded to no-data to reduce noise in the following analyses. We selected 44 European localities with different values of habitat suitability (See Table 1), and applied recursive partition trees (rpart library, \cite{rpartcitation}) using habitat suitability as response variable to group them in three different ways: 1) using the values of the environmental variables as predictors, to group the localities according to their environmental similarity; 2) using the local R-squared as predictors, to group the localities depending on the importance of the predictors; 3) using the local coefficients as predictors, to group together the localities with similar effect of the environmental variables.