5.1 Wind fetch

Wind fetch is a fundamental concept in the study of ocean systems and environments that is commonly employed outside archaeology to estimate the exposure of locations to wind-generated wave action (Laing 1998, for archaeological applications see Nitter, Elvestad & Selsing 2013; Nitter & Coolen 2018). Wind-generated wave action is a complex phenomenon that can depend on a wide range of parameters. However, the most commonly applied models for wave exposure are founded in so-called fetch-based indices. Fetch is defined as the unobstructed distance of sea across which the wind travels towards a given location. The longer the uninterrupted distance the wind can travel the surface of the ocean, the higher waves will tend to be. Wave prediction models based on fetch involve various degrees of complexity, starting from relatively simplistic models that only take account of the length of the fetch, as is done here. While more complex models allow for more accurate predictions that also move beyond estimation of mere relative exposure (Malhotra & Fonseca 2007; Bekkby et al. 2008; Sundblad et al. 2014), investigations into the predictive power of fetch-based models have found that those based solely on fetch length still perform reasonably well in predicting the presence and absence of biological indicators of wave exposure (Burrows, Harvey & Robb 2008; Hill et al. 2010).

Here, estimation of average fetch length was based on the procedures presented by Ekebom, Laihonen and Suominen (2003) and Tolvanen and Suominen (2005). With the sea-level raised just above each random point and the centre point for the site polygons, this involved first drawing a line every 7.5 degrees in circumference around each point. The lines were then clipped by the raised coastline in the study area and by the present day coastline in the larger region (Figure 3). The average length of the lines was then returned. The coastline has also changed considerably outside of the study region throughout the Mesolithic, which means that the average fetch of the more exposed locations has presumably been underestimated. However, as the average fetch lengths will largely be determined by the immediate surroundings of a location, this seemed like less of a problem for the purposes of mapping general locational patterns.

5.2 Visibility

Viewshed analysis was employed to investigate visibility patterns (Conolly & Lake 2006: 225–233; Gillings & Wheatley 2020). As my interest here were general visibility patterns, simple measures of total visible area from the analysed locations was deemed sufficient. The outer radii of the viewsheds was set to 5000 and 10,000 meters to represent views at a long and medium distance. These distances are somewhat arbitrary, but relatively common when mapping the size of views where a target has not been defined (e.g. Lake & Woodman 2000; Lopez-Romero 2008; Garcia 2013). The elevation of the observer point was set to 1.65 m above ground, and r.viewshed, the employed GRASS module, was set to take the curvature of the earth into account.

The sea level was raised to the lowest elevation of each site polygon, and the viewsheds computed from the centre points. For raising the sea level to the random point locations, a circular buffer with a radius of 11.9 meters was generated around each point. This equates to the median site size. Results were returned as the number of raster cells visible from each point.

5.3 Shoreline displacement

Wren, Costopoulos and Hawley (2020) recently investigated the relationship between settlement location of the Wemindji and rapid isostatic rebound in James Bay, Canada. They found that the sites in the region tend to be situated at locations where the topography will have resulted in a higher than expected degree of coastal stability within a 2 km radius of the sites, and a higher than expected variability in shoreline emergence within 20 km of the sites. This is taken to reflect a settlement strategy favouring site locations with stable local surroundings that simultaneously offer ecological variability within wider resource catchments. Variation in when areas emerged from the sea means that these will be at various stages of ecological succession, following the transition from oceanic, to littoral, to terrestrial environments.

A similar approach is taken here, with a few adjustments to the methodology of Wren, Costopoulos and Hawley (2020). First, the shoreline displacement curves for the study region are only reasonably reliable up to around 100 masl. The 32 EM sites situated above 100 masl were consequently excluded from this part of the analysis. Furthermore, the shoreline displacement in south-eastern Norway cannot be approximated by a low-order polynomial. The maximum radius explored around sites and non-sites was therefore set to 1 km, due to uncertainty and variable uplift rates towards the outer edges of the distribution of sites.

The emergence raster in Figure 4 was created by reclassifying the elevation raster based on the regional displacement curves. As the assumption here is that all sites have been situated on the shoreline, all areas lower than the lowest elevation of the site polygons and non-site buffers were excluded from analysis for each location. Standard deviation in year of emergence was then retrieved using buffers with radii of 50, 500 and 1000 meters around the point locations.

5.4 Island location

Several authors have pointed to the fact that a substantial number of sites seem to have been situated on small islands, and especially so in the earlier parts of the Mesolithic (e.g. Nyland 2012; Bjerck & Zangrando 2013; Bjerck 2017). Behavioural explanations of this assertion has followed a general emphasis on the use of watercrafts and resource exploitation in archipelago environments, but has also been explicitly related to extensive seal hunting in the earliest part of the Mesolithic (Bjerck 2017), as well as concentrations of marine productivity in the lee of smaller islands (Breivik 2014; Schmitt 2015). For this analysis a context dependent algorithm available in the provided Python script was devised to identify if point locations would have been situated on islands, assuming they were shore-bound. If yes, then the size of the island was found.

5.5 Sediments and infiltration ability

It has also been suggested that Mesolithic sites tend to be located relative to different soils and sediments (see Bergsvik 1995; Berg-Hansen 2009). This suggestion is based on the assumed desire of hunter-gatherers to situate sites on well-drained locations. The sediment data used here was retrieved from the Geological Survey of Norway and has a resolution of 1:50,0000. The sediment data has an infiltration score, relating to how well the sediments drain and clean wastewater. While these are in part related to the potential for biological breakdown of toxins, they are also determined by degree of permeability (Geological Survey of Norway 2015). Thus, while the original motivation for the classification does not perfectly align with my intended use, it does appear to capture the relevant properties of the sediments. The infiltration level for each location was found by taking the mode on the site polygons and the non-site buffers.

5.6 Aspect

Aspect may also have been important to locational choices in Mesolithic Norway. This is typically seen in relation to shelter (Berg-Hansen 2009: 47), a desire to face either land or the open sea (Darmark, Viken & Johannessen 2018), or to situate sites at locations that receive more sunlight throughout the day (Mikkelsen 1989: 58; Berg-Hansen 2009: 113). Here it was decided to focus on the possible relevance of aspect as related to sunlight, as the other proposed reasons for why aspect might be relevant should be captured by other variables. The aspect of each location, calculated as degrees from north, was therefore transformed to deviance from south by finding the absolute distance of each value from 180 degrees (Spencer & Bevan 2018).

6.1 Model building strategies

The first statistical method applied is logistic regression, which forms the backbone of site locational analysis in archaeology. The most common approach to model building with regression methods in archaeology could be said to fall in the category of stepwise minimisation, sometimes also called data-driven model tuning, as this often starts with analysis of the relationship between predictors and the response, independently of the multivariate model (Conolly & Lake 2006: 182–183). This becomes very problematic when variables are included or excluded from further analysis based on the statistical significance of this relationship, as has sometimes been the case (e.g. Stančič & Veljanovski 2000: 152–153; Duncan & Beckman 2000: 36). This is in effect a form of p-value driven, forward stepwise variable selection, and is hampered by severe issues (Harrell 2015: 71–72). First of all this approach would exclude a variable without considering whether it might be an important predictor once other variables are controlled for, or if it might alter the effect of other predictors in the model. If such effects are present, the variable should be included in the model irrespective of whether or not it is a good predictor of the independent variable (Kohler & Parker 1986: 416; James et al. 2013: 89). Secondly, this method has a tendency to inflate model performance. This follows from the fact that the final model specification is typically evaluated using standard methods that do not take the model tuning process into account. Not accounting for the total number of variables considered, including those that were discarded, breaks the assumptions underlying most standard statistical hypothesis testing and estimation methods, leading to a deflation of p-values and confidence intervals, and an inflation of R² values and coefficients (Harrell 2015).

However, the term stepwise selection is typically taken to refer to a form of automated model selection that has also seen application in archaeology (e.g. Bevan et al. 2013; Visentin & Carrer 2017; Spencer & Bevan 2018; Wachtel et al. 2018). Here, independent variables are excluded or included in the model at sequential steps that are stopped once a statistic that determines the relative success of the models does not improve. The Akaike and Bayesian information criteria (AIC and BIC, respectively) are the most frequently applied statistics for stepwise model selection in archaeology. Stepwise selection based on AIC and BIC are seen as an improvement over simple minimisation by use of p-values as a stopping criteria, as they take the total number of variables considered into account in the evaluation of the final model. The fundamental strategy of automated stepwise model selection has, however, received considerable criticism (e.g. Henderson & Denison 1989; Derksen & Keselman 1992; Chatfield 1995; Malek, Berger & Coburn 2007; Burnham, Anderson & Huyvaert 2011; Harrell 2015), and the implications of multiple testing and model tuning are more profound and difficult to take account of in a satisfactory way, leading Harrell (2015: 69) to conclude that ‘AIC is just a restatement of the p-value, and as such, doesn’t solve the severe problems with stepwise variable selection other than forcing us to use slightly more sensible α values.’

The most basic model building strategy associated with regression techniques is the forced entry, direct or simultaneous technique, where all independent variables are included from the start, with no consideration of the relative importance, or the order by which the variables should be included (e.g. Stoltzfus 2011). This is typically best suited when there is little reason to hold any a priori assumptions about the nature of these aspects, and is the approach taken here. The main drawback of this approach is the danger of overfitting and including irrelevant variables that only contribute an increase in standard errors.

Given the dramatic increases in computational power over the last decades, iteratively fitting hundreds or even thousands of models has now become tractable. Faced with this, automated model specification through stepwise selection could appear to be one way to handle this flood of information. Furthermore, regression models can offer insight into the complex net of relative and absolute variable importance, the relationship between these, and an estimate of the confidence with which these assertions can be made. But, as shown above, several authors have argued that this potential is lost with the implementation of stepwise model selection. To fully utilise their capacity arguably requires decisive and rigorous implementation of a well-thought-out, thorough and at times time-consuming modelling strategy (e.g. Harrell 2015). It could perhaps be argued that in the face of a fragmented and sparse material such as that frequently encountered in archaeology, a looser and more exploratory modelling procedure is warranted. Woodman and Woodward (2002) have argued the precise opposite. They contend that this uncertainty is precisely where careful modelling is most necessary and beneficial, as this can help both elucidate the nature of the fragmentation and guide an analysis that could otherwise be led astray.

Relevant to this is the delineation that Breiman (2001b) makes between what he terms a data modelling culture and an algorithmic modelling culture in statistics. The data modelling culture is based on an assumption that the data under study represents independent draws from a data-generating process where the value on a response variable follows a stochastic model of the form f (parameters, predictor variables, random noise). The algorithmic modelling culture, on the other hand, only makes the assumption that the data represents independent results of a complex and mostly unknowable underlying multivariate process. The goal is to identify an algorithm that best predicts the outcome of this underlying process. Good prediction indicates a better emulation of the underlying process than a solution that predicts poorly. A normal criticism of this focus is that prediction is not explanation, and that despite predictive success, the algorithms in use are often too impenetrable to contribute much in the way of explanation. Breiman (2001b: 210) states that this dichotomy of interpretability and explanation as opposed to prediction is contrived: ‘The goal is not interpretability, but accurate information’. While regression models can offer a highly interpretable output in the form of coefficients and confidence intervals, they are argued to be easily distorted by resting on assumptions concerning the nature of the underlying distributions that are virtually never met in the analysis of real world phenomena. Additionally, they appear to often be outperformed when it comes to prediction. While the output of machine learning techniques is less interpretable, it is certainly less opaque than the original process under study. Hence, it is argued, they can provide far more reliable and accurate insight, even if this insight holds less detail.

6.2 Random forests

Relatively high predictive power, ease of implementation and its non-parametric nature are among the elements that have led to the popularity of random forests (Hastie, Tibshirani & Friedman 2009: 587–604). Random forest is an ensemble method that leverages the increased predictive power gained from aggregating the results of a large number of de-correlated, high-variance, low-bias decision trees (Breiman 2001a). Decision trees for classification are based on identifying what variables best partition the data into classificatory groups at sequential splits (Baxter 2003: 116–118). At each split the variable to use is decided based on how well the candidate variables partitions the data into the correct groups. The degree of discriminatory success is also termed purity, and is often given by the Gini impurity score, which represents how often a random observation would be labelled incorrectly on the given split. While classification trees are highly interpretable and often perform well in classifying the original training data, they have a tendency towards high variance (Hastie, Tibshirani & Friedman 2017: 312).

To avoid overfitting while reducing the amount of variance, random forest works by combining bootstrap aggregation of the training data with a random selection of predictor variables to be used in the construction of each individual tree (Box 1). Thus, for each tree the training data is first randomly redrawn with replacement, which results in around 2/3 of the original data being evaluated in the training of the entire forest—a process also known as bagging. The remaining 1/3 is known as the out-of-bag (OOB) sample, and can be used for subsequent model validation and variable importance estimates. Following each resampling, a classification tree is grown. For each node in the tree, a random subset of predictors is drawn, and the predictor returning the lowest Gini impurity score is selected at that node split. The number of randomly selected predictors evaluated at each node is recommended to start at $\sqrt{p}$M4 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \sqrt p \] \end{document} for classification trees, where p is the number of independent variables. The number of variables to be evaluated at each split, the parameter mtry, is then explored during model tuning. As there is no danger of overfitting associated with generating too many trees, the arbitrarily large number of 5000 trees were generated for each forest here. While this came at a computational cost, it meant that only the mtry parameter required tuning.

An important output of random forests is the variable importance estimates that rank the variables based on their ability to correctly classify the input data. Following Parr et al. (2018), variable importance was here found by means of the permutation importance method. This is based on first running the OOB-sample through the final forest to establish the baseline accuracy. Then the column of values for an independent variable is randomly shuffled, and the OOB-sample is run through again. The difference in accuracy between the two runs then gives the variable importance measure. It is crucial to note that these can only be used to evaluate the relative importance of variables, and do not provide estimates of effect size.

6.3 Model validation

Internal model validation was done by means of bootstrap resampling for the logistic regression models (Kvamme 1988; Verhagen 2009). The purpose of applying this method is to replace distributional assumptions and asymptotic results with computationally derived confidence intervals (Fox 2008: 605), and to avoid both fitting and evaluating model performance using the exact same data (Hastie, Tibshirani & Friedman 2017: 249–254). This entailed iteratively resampling the original data and fitting each model 9999 times. For each iteration the coefficients were retrieved, in addition to the Brier score and Area under the receiver operating curve (AUC). The Brier score is a proper accuracy score, meaning that it summarises calibration, the difference between the predicted probability and the observed probability of an outcome, as well as model resolution, the spread of the predictive distribution (Rufibach 2010). The measure ranges from 0 to 1, where 0 represents perfect accuracy. The score is not especially meaningful on its own, however, and is best suited for comparative purposes. The AUC score has a more intuitive interpretation than the Brier score, but only considers the discriminatory abilities of the model (Hosmer, Lemeshow & Sturdivant 2013: 173–182). The AUC ranges from 0 to 1, where 1 represents perfect classification, and 0.5 represents no improvement over random classification. The accuracy scores are reported as the mean result across bootstrap samples. The coefficients were evaluated by means of 95% bootstrap percentile intervals, constituted by the range between the 2.5th and 97.5th percentiles of coefficient values generated across the bootstrap samples (Fox 2008: 595).

Variable importance estimates and model performance measures for the random forests were found using nested k-fold cross-validation instead of the bootstrap (Box 2, and below for the imputation of missing values). The reason for this is that while the logistic regression models were fit using the simultaneous technique, a model tuning step was included for the random forests to identify the optimal number of random variables to be evaluated at each split in the classification trees (the parameter mtry). Single bootstrap resampling results in some overlap between the original data and the resampled observations, meaning tuning and validating with this technique would likely lead to overfitting and overestimation of model performance (Hastie, Tibshirani & Friedman 2017: 250). Basic cross-validation involves randomly splitting the data into k folds, training the model using all but one fold, and then validating the result against the retained fold (Verhagen 2009: 92; Hastie, Tibshirani & Friedman 2017: 241–249). This is then repeated until each fold has been used as the validation set. The nested version of cross-validation entails that for each of these validation steps the folds not retained are split into an additional k folds. Here a complete sweep of possible mtry settings was done for each inner cycle, and the mtry achieving the best performance by means of AUC was returned and used for the current outer fold. The final AUC and Brier scores, as well as variable importance measures are reported as averaged across the outer validation cycle.

6.4 Model parameters and pre-processing

Below follows a presentation of predefined model parameters and the pre-processing steps that were taken. With the exception of imputation, the response variable was kept out of this process to maintain inferential validity.

Apart from aspect deviating from south, all continuous variables were transformed for the logistic regression models by taking the natural logarithm due to extreme distributions. Continuous variables were then normalised to take on values between 0 and 1 to better facilitate a comparison of the impact of each variable, although at the cost of distorting the interpretability of effect sizes (cf. Spencer & Bevan 2018). In the random forests, the variables were included using raw values. For a second set of models directly comparing the location of sites between the different phases, as opposed to sites and non-sites, all variables with the exception of the binary island/mainland variable were transformed by finding the percentile rank of each site value as compared to the values in the corresponding hull sample. In an attempt to achieve more sensible diachronic comparisons, this adjusted the values by variation in the surrounding landscape as captured by the hull samples for each phase.

There are a host of different methods for identifying and dealing with collinearity (e.g. Dormann et al. 2012; Tomaschek, Hendrix & Baayen 2018). The approach taken here was a relatively straightforward one, where Pearson’s r and Spearman’s ρ was found for each predictor variable pair for each chronological phase. In the case of problematic levels of correlation between predictor variables across all phases (|r| or |ρ| > 0.8), all but one of the correlating variables were removed from analysis. The choice of what variable to retain was based on what variable appeared to best summarise the collinear group in substantive terms (Table 2).

Following the lack of previous quantitative studies on which to base this analysis, linearity with the log-odds of the response could not be assumed for any of the predictors in the logistic regression models. To allow for non-linearity, pre-specified natural splines with default quantile knots could have been included (Harrell 2015: 26–28). However, with 104 sites for the EM representing the smallest effective sample size, the inclusion of more than 10 model parameters would likely lead to overfitting, following the most lenient m/10 rule provided by Harrell (2015: 72–73), where m is effective sample size. Consequently, as there were no grounds on which to prioritise among variables to include with splines, none were used. Similarly, due to the large number of possible two-way interaction effects, no interaction terms were included in the models as there were no grounds on which to prioritise among these (cf. Harrell 2015: 36–38).

In total, shoreline emergence was missing for 32 sites and 725 non-sites, while infiltration level was given as unclassified for one site and 138 non-sites. Imputation of these missing values was deemed most appropriate. Unless the number of cases that have missing data is very low, or the values are missing completely at random, which is rarely the case, the alternative of deleting observations can lead to dramatic loss of information and bias in all model estimates (Nakagawa & Freckleton 2008). Imputation for the logistic regression models was done using predictive mean matching (van Buuren & Groothuis-Oudshoorn 2011). This is based on performing ordinary least squares regression on the column with missing values, using the other variables in the data as predictors. From the resulting coefficients a random draw is then made, and these coefficients are used to predict the entire column containing missing values, including those values that are not missing. Each missing value is then given the originally observed value of one of its predicted closest neighbours by a random draw, typically among the five closest cases. This process is then ideally repeated multiple times, and subsequent analysis is based on a pooling of the results from running each imputed dataset through the full analysis. As the imputations done for the logistic regression were to be nested in a bootstrap, a single imputation was performed per bootstrap sample (Brand et al. 2018).

Random forest provides its own method for imputing missing values based on its proximity and classification capabilities (Breiman 2003). This works by first setting the missing values to the average of the column and then growing and running the dataset through a random forest. For each classification tree in the forest, each observation that ends up in the same terminal node as the observation with the missing value is counted as similar. The value to be imputed is then given as the average across the proximal observations, weighted by the number of times they ended up in the same terminal node. This entire process is then repeated using the new values instead of the column average. Here five sequential random forests were grown, with each forest growing 5000 trees. This process was then repeated five times, creating five imputed datasets that each were run through the estimation and validation process described above.

7.1 Settlement patterns of Mesolithic hunter-fisher-gatherers

A shift from a resource base focused on the hunting of sea mammals towards a wider subsistence spectrum that includes the regular exploitation of ungulates, birds, shellfish and especially an intensification of fishing is a gradual development commonly believed to characterise the Norwegian Stone Age (Breivik 2014; Bjerck et al. 2016; Ritchie, Hufthammer & Bergsvik 2016; Jørgensen 2020; Mjærum & Mansrud 2020). Furthermore, a transition to a more diverse marine economy centred around fishing has been argued to have major implications for the mobility patterns and settlement strategies of coastal hunter-gatherers, where such a shift could indicate a decrease in mobility (see above and Boethius 2017; Boethius & Ahlström 2018). The location of MM and LM sites with respect to variation in shoreline emergence, as identified in this study, could therefore be consistent with economic diversification and a decrease in mobility. However, the negligible performance measures and effect sizes in the comparison between phases given in Figures 12, 13, 14 do not allow for the conclusion that the settlement patterns underwent any consequential changes between the three considered phases. In addition, the models comparing sites and non-sites indicate that location relative to exposure has for the most part been more important than the shoreline emergence variables. Furthermore, the effect sizes for the two exposure variables are comparable across all phases, and in the comparison between phases these variables have a relatively steady impact around zero, indicating that site location with respect to exposure has been similar in all phases. Thus, if degree of site exposure is taken to be a reflection of the economic basis for settlement, as has been proposed (e.g. Bjerck 2009; Breivik 2014), the patterns discerned here would indicate a stable subsistence base throughout the Mesolithic. An alternative explanation to the apparent stability in settlement patterns could of course also be that these might simply not be sensitive to changes in the subsistence patterns of coastal hunter-fisher-gatherers in the archipelago environment of south-eastern Norway, and that this setting demanded that settlement patterns take on a distinct form, irrespective of variation along other social and cultural dimensions. One possible reason for this could perhaps follow from a dependency on boating that this landscape is likely to have represented (e.g. Bjerck 2008; Glørstad 2013). This might, in turn, have constrained the variation space available for locational strategies. It is, however, important to underscore that the limited performance of the models do mean that while no conclusion of clear diachronic difference can be drawn, this does not mean that the settlement patterns were necessarily the same—a considerable amount of variation does remain unaccounted for.