South Arabia’s prehistoric monument landscape shows social resilience to climate change

PLOS One | https://doi.org/10.1371/journal.pone.0323544 May 28, 2025 8 / 24 spring sites [80]. The distributions of travertines attest to former springs, and spring output today is lower than in the past. Springs that persist today in the Nejd were likely also flowing in the past and were therefore particularly reliable water sources around which we expect human activity to be high. Vegetation was measured by calculating Normalized Difference Vegetation Index (NDVI) from all available imagery in the Landsat 8 surface reflectance imagery collection in Google Earth Engine after removal of very cloudy images (>30%) and cloud-masking. A map layer of the maximum NDVI values was used for analysis rather than NDVI because values are generally very low and highly variable in arid environments. Fog density quantifies the average distribution of monsoon fog in the study area. The layer was calculated on a per cell basis as the mean of the near-infrared fog reflectance values of 119 Landsat 5 TM scenes, 17 Landsat 7 ETM+ scenes and 121 Landsat 8 OLI TIRS scenes [81]. Maximum NDVI and fog density correlate with a major bio-climatic zonation in Dhofar, which persists due to the inter- play between topography and monsoonal advection fog [81,82]. The plateau and the escarpment would always have been moister than the interior and form a natural axis against which to measure regional distributions in the arid interior. We also include eight topographic variables in our models because topographic factors such as aspect, visibility and/or acces-

sibility may have influenced more localized monument placement [59,67]. We evaluated variables similar to those tested elsewhere, including in spatial modelling of ancient funerary monuments in eastern Sudan [2]. Bootstrap aggregating. Bootstrap aggregating (bagging) is a machine learning method that averages the results of multiple classification trees, each trained on different random (bootstrap) samples. Combining the results of multiple models helps to improve accuracy and reduce overfitting, compared to using a single tree. Combining the results of multiple iterations of models is common and helps improve accuracy by balancing out each model’s strengths and weaknesses, leading to better, more reliable predictions. Classification trees split data into smaller subgroups through recursive partitioning, based on true or false answers about the values of predictors. Each split is based on a single variable, and a rule that maximizes the number of similar values (class purity) within each of the two resulting subgroups [83]. For example, a classification tree might split animals by weight, classifying those heavier than 10 kg as mammals and those lighter as birds. Classification trees are well-suited to geospatial datasets where nonlinear responses and predictors may interact in unknown ways [84]. In bootstrap aggregating, classification trees are grown until all the observations have been classified, unlike single trees where a complexity parameter limits tree growth. Bootstrap aggregating was

performed in the R ‘rpart’ package [85,86] which implements methodologies of Breiman et al. [83]. We chose to generate a high number of trees (500) as the results of single trees differed greatly. Model accuracy is assessed using the out-of- bag misclassification error. Variable importance was calculated as the sum of reductions in mean-squared error each time a variable was used as a primary or surrogate split and summed across all trees. For each variable, partial dependence plots were generated to visualize the relationships between the environmental variables and the monument locations [87]. (S1 Fig) Environmental data analysis and regularized multinomial logistic regression To examine the power of the environmental variables for predicting monument types and to understand the direction and magnitude of these relationships, we performed regularized multinomial logistic regression. This method was preferred over machine learning techniques such as bootstrap aggregating, as it is more robust to the unequal numbers of observa- tions of different monument types. It also incorporates variable selection and provides highly interpretable results. Reg- ularized multinomial logistic regression was performed in the ‘glmnet’ and ‘glmnetUtils’ packages in R [85,88]. A range of regularization techniques can be applied, from ridge regression [89], to elastic net [90], to Least Absolute Shrinkage and Selection Operator (LASSO) [91]. Parameter α determines the type of regularization. When α = 0, ridge regression shrinks

the coefficients of correlated variables toward each other but retains them all. When α = 1, LASSO shrinks coefficients, producing a sparse model by filtering out less important variables. When α = > 0 or < 1, elastic net regression balances