Journal of Archaeological Science: Reports 11 (2017) 59–68 Contents lists available at ScienceDirect Journal of Archaeological Science: Reports journal homepage: www.elsevier.com/locate/jasrep A comprehensive test of the Locally-Adaptive Model of Archaeological Potential (LAMAP) W. Christopher Carleton a,⁎, Kong F. Cheong b, Dan Savage c, Jack Barry c, James Conolly c, Gyles Iannone c a b c Human Evolutionary Studies Program and Department of Archaeology, Simon Fraser University, 8888 University Drive, Burnaby, British Columbia, V5A 1S6, Canada Department of Anthropology, American University, 4400 Massachusetts Ave NW, Hamilton Hall, WA, District of Columbia 20016-8003, United States Department of Anthropology, Trent University, 2140 East Bank Drive, Life & Health Sciences, DNA Building, Module “C”, Peterborough, Ontario K9J 7B8, Canada a r t i c l e i n f o Article history: Received 18 June 2016 Received in revised form 9 November 2016 Accepted 15 November 2016 Available online xxxx Keywords: Classic Maya GIS Predictive modeling LiDAR Survey Human-landscape interaction a b s t r a c t Extensive archaeological surveys are critical for understanding past human-landscape interaction, but they are frequently impeded by access difficulties, rugged terrain, or obscurant vegetation. These challenges can make extensive surveys prohibitively costly and time-consuming. Consequently, many archaeologists are interested in predictive techniques—i.e., methods that can estimate the potential for a given region to contain archaeological remains. Predictive techniques can reduce the costs of extensive surveys by allowing archaeologists to focus on the regions with the greatest archaeological potential. A few years ago, our research team developed a new technique called the Locally-Adaptive Model of Archaeological Potential (LAMAP) and used it to enhance our understanding of the relationship between the Classic Maya centre of Minanha, its surrounding landscape, and nearby Maya centres (Carleton et al. 2012). However, when we introduced the method its efficacy had yet to be comprehensively tested. Recently, we tested its efficacy using a combination of ground-truth survey and remote sensing of Classic Maya sites in west-central Belize. The test involved identifying previously unrecorded archaeological resources and comparing their locations to the LAMAP prediction and to a random model that acted as a null hypothesis. Our results indicate that the model performs very well. The high-potential areas of the study region contained three times more archaeological sites than low potential areas, a statistically significant result compared to our null model. Our findings indicate that the LAMAP is a useful new archaeological prediction tool and, as a corollary, that the hypothesis of human land-use behaviour underpinning it might accurately reflect the behaviour of the Classic Maya. © 2016 Elsevier Ltd. All rights reserved. 1. Introduction Human-landscape interaction is a foundational element of human ecology. It has defined and constrained where people could live, how they could feed themselves, how they could communicate or move about, and the natural resources that they had access to, affecting our biological and cultural evolution. Thus, finding patterns in past humanlandscape interaction can illuminate important aspects of human history, including evolutionary processes, migratory events, economic and political structures, and even past human cognition (e.g., Alexakis et al., 2011; Bevan and Conolly, 2002; Conolly and Lake, 2006; Kanter, 2008; Llobera, 2003; Rua, 2009; Swanson, 2003; Verhagen and Whitley, 2011). To see these patterns we often have to use extensive archaeological surveys because they provide the necessary data at a large ⁎ Corresponding author. E-mail addresses: w.ccarleton@gmail.com (W.C. Carleton), kong.cheong@american.edu (K.F. Cheong), danielsavage@trentu.ca (D. Savage), jack.barry29@gmail.com (J. Barry), jamesconolly@trentu.ca (J. Conolly), giannone@trentu.ca (G. Iannone). http://dx.doi.org/10.1016/j.jasrep.2016.11.027 2352-409X/© 2016 Elsevier Ltd. All rights reserved. enough scale to observe correlations between human activity and landscape features or traits. They can reveal, for example, patterns in the proximity of settlements to various natural resources or differences among settlements in terms of ecological conditions, from which we can infer past landscape use. Large-scale surveys can be challenging, however, because of difficult access, rugged terrain, dense vegetation, and sediment accumulation. These challenges can make traditional pedestrian surveys prohibitively time consuming and expensive. And since archaeological projects are commonly constrained by limited resources, complete surveys of large areas are not always possible limiting our view of large-scale patterns in past landscape use. The challenges involved in large-scale landscape survey can be partially overcome with predictive modeling. Predictive models are tools for estimating the potential of a given parcel of land to contain archaeological material on the basis of a set of empirically observable relationships between site locations and a set of potential predictor variables (Verhagen and Whitley, 2011). In the last few decades, several predictive modeling techniques have been published in the academic literature, but rarely have they been rigorously tested after their initial publication. Here we report the results of a robust test of a new 60 W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 predictive modeling technique called the Locally-Adaptive Model of Archaeological Potential (LAMAP; Carleton et al., 2012). We developed the LAMAP method to study past landscape use around the Classic Maya centre of Minanha located on the northern Vaca Plateau in west-central Belize. The region around Minanha is characterized by exactly the type of terrain that can impede large-scale surveys. The plateau is remote, difficult to access, extremely rugged, and covered with dense jungle—an eminently suitable place to develop and test predictive modeling techniques. 2. Background Now a common tool in archaeology, predictive modeling has advanced considerably in the last few decades (McCoy and Ladefoged, 2009; Verhagen, 2007; Verhagen and Whitley, 2011). The advances have been facilitated by the widespread adoption of Geographic Information Systems (GIS) among archaeologists and the increasing availability of digital topographic and environmental maps, high-resolution satellite imagery, open-source software, and powerful computers. These tools have been used to create predictive models for regions all over the world (e.g., Demján and Dreslerová, 2016; Dorshow, 2012; Graves, 2011; Howey, 2015; Oonk and Spijker, 2015). These predictions constitute an expectation about where unrecorded sites are most likely to be found—essentially, a hypothesis about past human-landscape interaction. Generally, however, archaeological predictive modeling is based on ‘presence-only’ data as true site absences are unknown, and thus predicting the true probability of site occurrence is impossible (Verhagen, 2007). Rather, they establish a prediction of relative suitability of a given parcel of land for past occupation. Often, the models are then used to generalize about past human landscape use or support heritage resource management and planning in the government and private sectors. The majority of archaeological predictive modeling techniques share common characteristics despite differing in their particulars. There are several established modeling techniques and the newer ones take advantage of modern Bayesian statistics, high-performance computing, and computer learning algorithms that are also used in applications like artificial intelligence (e.g., Fernandes et al., 2011; Hitchings et al., 2016; Oonk and Spijker, 2015). Most of these approaches involve regression models for estimating relationships between known-site locations and a set of predictor variables such as slope or distance to water. The regression equations are then used to calculate the probability that a given location has landscape traits that make it likely to have been favorable for past people as an activity area or settlement— i.e., a measure of its potential to host a site, rather than what is often and incorrectly described as “probability of site location”. These probabilities are generally derived for every cell in a raster image of a given region and classified into categories of potential like “high”, “medium”, or “low”. The efficacy of the model is usually assessed with a measure such as Kvamme's Gain Statistic, which is a ratio of the number of known sites located in highpotential cells to the percentage of the region classified as high-potential (Kvamme, 1988). Often the models perform well in the sense that they can “predict” the locations of the sites in a known-site database and some sample of that database that was withheld to validate the model. However, many authors have leveled important criticisms against the existing approaches (e.g., Ebert, 2004; Mehrer and Wescott, 2006; Wheatley and Gillings, 2002). Verhagen (2007) summarizes the various criticisms, condensing them into several areas that he argued need improvement in predictive modeling. First, he argues that the archaeological data used in predictive modeling is often too incomplete to provide a meaningful picture of past landscape use, leading to biased models. So, according to Verhagen, we ought to build up a more complete dataset in order to create better predictions. Second, he stresses that the spatiotemporal resolution of the archaeological and environmental data is not high enough to make good predictive models. Third, Verhagen states that we need to improve our statistical tool-box because our techniques are out of date and, consequently, so are the existing predictive models. Fourth, he underscores that most predictive models do not include sociocultural information, leading to biases in favour of environmentally deterministic views of past human–landscape interaction. And, lastly, he suggests that we ought to test our models more rigorously and more often. While we agree with Verhagen's assessment of the problems with predictive modeling in general, we think that there are other, more specific reasons that the existing approaches are problematic, some of which were explored by the Archeomedes project (e.g., Durand-Dastès et al., 1998; Van der Leeuw et al., 2003). The first problem is that they are overly focused on point data—i.e., treating archaeological sites like points instead of areas. This undue attention to point data reduces site locations to one-dimensional phenomena rather than local regions that encompass a settlement area of potentially diverse characteristics. And the specific characteristics of a given region that made it attractive to human settlement may be difficult to describe by summary statistics, like the mean or range. It may instead be the bior multi-modal features of the landscape variables that are attractive—e.g., multiple modes in elevation or aspect. A second problem is that predictive models over generalize. They typically apply the same predictive equation to every area within a given region. So, they cannot account for the potential geographic variation in site location preferences without creating separate models. For example, it is unrealistic to assume that characteristics favorable to settlement in a lowland wet area would be similar to those in a highland arid area—i.e., human-landscape interaction is often spatially heterogeneous. In ecologically complex and heterogeneous regions, the models are at risk of missing variation in land use strategies, conflating site types and giving a false impression of past landscape use. The same criticism applies to temporal differences. The third problem, to which we alluded earlier, is that archaeological data generally consist of presence-only observations, but many of the popular techniques are founded on having a more complete dataset. Archaeologists have to know where sites are not located along with knowing where they are located in order to properly use many existing methods. And for the same reasons that large-scale systematic surveys are challenging—namely, costs and access impediments—intensive random sampling to identify locations that are truly absent of archaeological remains is often impossible. A popular solution has often been to use pseudo-absence data (e.g., Tallavaara and Seppä, 2012), like sets of random locations, but the approach is problematic (Barbet-Massin et al., 2012; Chefaoui and Lobo, 2008). It technically violates the statistical assumptions of the established methods because the status of the archaeological record at the random locations is actually unknown, not negative. This problem is not unique to archaeology since ecologists and biologists interested in species distributions modeling also often have to deal with presence-only data (Hastie and Fithian, 2013), which prompted the development of GARP and MaxEnt models in those fields (Phillips et al., 2006). Although archaeologists have explored these new methods (e.g., Conolly et al., 2012; William et al., 2006), they arguably suffer from the over-generalization problem because they cannot account for behavioural non-stationarity across space, and they still depend on simulated absence data. Thus, the lack of real absence data continues to be a challenge for archaeological predictive models, one that undermines all of the methods we are aware of. The final problem is that the established predictive modeling techniques are atheoretical (Verhagen, 2007). Most of the methods used for predictive modeling have been adapted from other sciences, like geology and conservation biology, and they are based purely on associations between presence/absence locations and potential environmental covariates. Verhagen and Whitley (2011) describe these models as “correlative”, “inductive”, or “data-driven”. These approaches make purely empirical claims with no attempt to model the underlying process. As such, the models do not account for the W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 complexity of human decision-making and, consequently, cannot be used to better understand past human land-use behaviour except in a correlative way. Although there is a large and robust body of theory and empirical data concerning the relationship between global environmental characteristics and human settlement patterns (e.g., Binford, 1982, 1980; Bird and O'Connell, 2006; Smith and Winterhalder, 1992), at regional scales the link between expectations of human land use and resulting patterns of archaeological site locations are typically uncritically developed. To overcome these problems, we developed a predictive method called the Locally-Adaptive Model of Archaeological Potential (LAMAP; Carleton et al., 2012). The LAMAP technique addresses each of the problems. Firstly, it uses distributions of landscape variables, instead of relying on point estimates or ranges. Secondly, it can adapt to localized landscape variation in known archaeological site locations and variation in site types, including temporal variations, adjusting its predictions accordingly. Thirdly, it does not involve knowing a priori where archaeological resources are not located, making it suitable for the majority of archaeological work. And, lastly, it makes a specific theoretical prediction about past human land-use decisions. The LAMAP technique can be thought of in terms of the theory that underpins it and the calculations used to operationalize that theory. The LAMAP theory assumes that people have certain notions about how to use their landscapes. Those notions are the product of biological and cultural evolution—a long process of trial and error by preceding generations who learned how best to use their landscape for whatever purpose. The notions serve as archetypes of land-use, blueprints for where to perform certain activities, including everything from shortterm camping to establishing permanent settlements or engaging in ritual. In this way, the LAMAP theory has an element of “group agency”, a term used by Verhagen and Whitley (2011) to describe large, grouplevel patterns of behaviour and collective norms. People pass the archetypes on to their children and peers through social learning and preserve them in cultural memory. The archetypes can change over time depending on the needs of the people, their interactions with a dynamic environment, or their interactions with other groups. According to this theory, each time a person goes out into the landscape to engage in some activity, they decide where to perform the activity by comparing their available options to their mental archetype for the intended activity, selecting from the candidate locations the one that most closely matches the relevant archetype. The LAMAP theory incorporates variation in the landscape by assuming that the person making the land-use decision is going to select the best of the available options using archetypes from nearby areas preferentially. The preference makes sense since nearby locations are more similar to the candidate locations than distant locations would be. Thus, the nearby sites are more likely to be well suited to the terrain in the candidate locations than more distant sites would be, making nearby sites better blueprints. Ideally, operationalizing the LAMAP theory would involve comparing candidate locations directly to the set of archetypes, but the archetypes themselves are not directly accessible. They are notions in the minds of past people. So, instead, we assume that known site locations are realizations of the archetypes that, taken together, indicate what the past archetypes might have been. Calculating a LAMAP prediction involves several steps (Carleton et al., 2012). To begin, a predictive region is selected and a training dataset of known site locations is gathered. The prediction region is the area for which a LAMAP prediction will be calculated and the known site locations can be from anywhere. However, the model will obviously reflect archaeological potential more accurately if the training sites come from somewhere within the prediction region. Then, the potential predictor variables are selected—e.g., slope, elevation, distance to fresh water, distance to high quality soils, and so on. The values of those variables within spatial catchments around each of the known sites are extracted using a GIS, creating a multivariate distribution of landscape variables for each known-site location. The catchments allow us to base our 61 predictions on areas instead of point locations, an idea considered by other scholars as well (e.g., Verhagen et al., 2013). They can be of any size and shape as long as they are the same size and shape for every site—the shape and size is a variable that can be altered and examined.Next, every cell in an output raster map of the predictive region is compared to each known site to derive a list of probabilities. The probabilities describe the likelihood that a given point in the prediction region comes from a parent population similar to the multivariate distribution of landscape variables in the catchments around each known site. To estimate these probabilities, we used joint empirical cumulative distribution functions calculated in R, though other approaches might work as well (R Core Team, 2016). To illustrate the process, imagine a single location of interest (LOI) compared to a single known site on the basis of elevation. The known site has a distribution of elevation values, derived from the cells in a raster map that represent elevation in a catchment around the site. After rounding the values to a sensible level of accuracy—e.g., single meters of elevation, as we did originally—we then calculate the frequency with which each elevation value occurs in the catchment. These frequencies represent the probability of finding a spot in the area around the known site that has a given elevation. However, the probability of finding an exact elevation value is usually exceedingly low. So, instead, we use cumulative distribution functions (CDFs) and intervals around the LOI value. The CDF describes the probability of finding a given elevation or lower for the area around the site. It is used to determine how likely we are to find an elevation equal to or lower than the elevation at the LOI. We then select a reasonable interval around the elevation value at the LOI to calculate an upper and lower boundary, which provides us with the probability of finding that elevation plus or minus the size of the interval. Suppose that the LOI has an elevation of 105 m above sea level. Using an interval, say 100–110 m, we can check the CDF of elevation values around the known site to estimate the probability of finding an elevation of 110 m or less and the probability of finding an elevation of 105 m or less. The difference between the two gives us the probability of finding an elevation of 105 ± 5 m in the catchment around the site. This probability serves as a proxy for the similarity between the LOI and the known site in terms of elevation. To improve the proxy, the calculation is extended to include multiple variables by using joint-CDFs. The joint-CDF describes the probability of finding, for example, a spot that has an elevation of 105 m and that has a slope of 30% and that sits 500 m from the nearest source of fresh water, and so on. In the multivariate case, we would compare a vector—a list—of probabilities where each element of the vector is a distribution of values for a given landscape variable in a catchment around a known site. Then, we calculate the probability of finding a spot in that catchment with the same multi-dimensional values as the LOI, plus or minus a reasonable margin. The final value is the probability of finding a location in the catchment around the known site with landscape variables similar to the ones at the LOI. And by this we mean jointly similar. So, for a given LOI to have a high probability of being found in a given catchment, all of its landscape variables have to be similar to all of the variables at numerous locations within the catchment. This probability serves as a proxy for the likelihood that the LOI possesses whatever landscape traits made the known site location desirable. This process is then repeated for every site in the known-site database. The result is a vector of probabilities expressing the likelihood that a spot like a given LOI could be found in the catchments around each of the known sites. The calculations are repeated for every raster cell in a given prediction region, so each cell has a vector of probabilities with a length equal to the number of known sites. In the last step, the probabilities are combined to produce a single probability for each LOI—i.e., for every raster cell in a given prediction region. As part of this step, the probabilities are weighted by the distance between the LOI and the corresponding known site, giving more weight to nearby known sites. We used a smooth exponential decay 62 W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 function in the original LAMAP study, but others might be appropriate. This weighting is, in part, what makes the LAMAP prediction “locally adaptive” because nearby sites are more likely to represent mental archetypes with relevance to the local area. Lastly, the probability of the union of the vector elements is calculated. The probability of the union is the probability that the LOI could feasibly be found in the catchments around any of the known sites. To calculate it, we have to use the inclusion–exclusion principle (Andrica and Andreescu, 2009), which allows us to account for the intersections in the probability spaces. The final result is a value between 0 and 1 for each cell in a raster map of a given prediction region. This value indicates the archaeological potential of a given cell with higher values indicating higher potential. The estimate, however, should always be considered relative since the question being answered is what is the archaeological potential of a given location compared to many others. Thus, the estimates should then be scaled to the highest and lowest LAMAP values in a given area because “high-potential” is always referring to the potential of one location relative to a surrounding space, not an absolute benchmark. When we developed the LAMAP method (Carleton et al., 2012), we used a case study to demonstrate the technique. The case study involved an area around a Classic Maya urban centre called Minanha, located in west-central Belize (Fig. 1). The site is situated atop the northern extent of the Vaca Plateau, an elevated region extending from the Maya Mountains to the Belize River Valley. Minanha was a small city-state that emerged as a significant micro-regional power in 8th century CE (Iannone, 2010, 2005; Iannone et al., 2014a, 2014b; Schwake and Iannone, 2010). The Minanha research project determined that the centre was peripheral to larger nearby urban centres, like Caracol, for most of its history. In some ways it was an archetypal Classic Maya urban centre characterized by a civic-ceremonial core occupied by elites and a surrounding settlement supported primarily by maize agriculture. One of the research aims of the Minanha project was to better understand the relationship between Minanha, its surrounding landscape, and other nearby centres. This objective led to the development of the LAMAP method. We created a high-resolution LAMAP model with 30 × 30 meter pixels for a 280 km2 region around Minanha using the locations of 69 Classic Maya civic-ceremonial centres from the Electronic Atlas of Maya Sites (http://mayagis.smv.org/). The site database was first cleaned by removing obvious duplicates and correcting the locations of well-known centres. Any remaining locational errors or duplicates would only serve to make our present analysis more conservative since they would likely create predictive errors, decreasing the utility of the LAMAP predictions. The 69 sites that remained were mostly large urban centres with typical features like pyramids, plazas, and domestic buildings. We based the LAMAP model on several terrain variables, including slope, ruggedness, distance to water, and an index of soil arability. A portion of the known sites were withheld from the model initially and used to validate it. The LAMAP model appeared to have predictive utility using standard benchmarks like Kvamme's Gain. However, the validation was not a comprehensive test, leaving open questions about the LAMAP's real utility for resource management and for improving our understanding of past human land-use behaviour. To test the LAMAP method more comprehensively, we gathered new survey data from the Minanha study region. The survey data allowed us to test the hypothesis that the LAMAP method is useful for estimating archaeological potential, meaning that it should accurately reflect archaeological potential in the Minanha study area. If the hypothesis were supported, we reasoned, a corollary of it would be that the LAMAP theory reflects human land-use behaviour—more specifically that people compare candidate areas on the landscape to one or more sociocultural archetypes in order to decide where to conduct a given activity. 3. Testing the LAMAP method Our analysis of the LAMAP method proceeded in two stages. First, we conducted two surveys of the Minanha study area. The first was a pedestrian field survey conducted by Barry and Cheong (2012). This survey was carried out using a sampling grid with a 500 m x 500 m resolution (Fig. 2). The grid was set up over the LAMAP prediction area and a random sample of 100 survey blocks was selected for pedestrian survey. The blocks were randomly sampled from the original LAMAP prediction surface, stratified by LAMAP value quintiles calculated in GRASS, resulting in 20 blocks per quintile placed randomly in the study area. However, the field conditions made it impossible to fully survey each of the 100 blocks. The density of vegetation, the ruggedness of the terrain, and the lack of direct access with a truck on service and logging roads made a complete survey impossible, underscoring the need for predictive modeling efforts. Fig. 1. Map showing study area. The dark grey region corresponds to the LAMAP model. W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 Fig. 2. Map showing the LAMAP model. Lighter areas have higher potential. Survey grids visited by Barry and Cheong are represented to scale by white squares. By the end, the survey team visited a total of 14 blocks, covering roughly 350 ha (see Barry and Cheong, 2012 for more details). The blocks were visited by targeting coordinates within the square with a handheld GPS, getting as close as possible in a truck and then walking to the location. The blocks targeted were selected by the survey team while in field with the aim of maximizing the number of blocks visited given access and time constraints. The 14 blocks the team was ultimately able to reach contained dozens of features, including 34 patio groups, hundreds of agricultural terraces, and three new minor centres. Two of the newly identified centres—named Kolchikiin and Ixkuk by the survey team—were on terrain deemed high potential by the LAMAP model, corresponding to classes 4 and 5 in Table 1 below. The third centre was located on terrain corresponding to LAMAP classes 3 or lower. Unfortunately, too few blocks were surveyed to be statistically informative. So, we excluded the ground survey results from the quantitative analyses that follow. Nevertheless, we were able to use these findings to supplement the next survey. The next survey was conducted in the lab using high-resolution LiDAR data of the study area collected in 2013. A portmanteau of the words Light and the acronym RAdio Detection And Ranging, LiDAR is one component of a larger technological system that enables the creation of high-resolution surface models. Over the last decade, LiDAR has Table 1 LAMAP classes and quantile ranges. LAMAP class Quantile low bound Quantile high bound 1 2 3 4 5 0 0.014049 0.056869 0.162760 0.347864 0.014049 0.056868 0.162759 0.347863 1 63 become an important part of archaeological research, with many recent examples of its successful application (e.g., Bewley et al., 2005; Challis, 2006; Chase et al., 2012, 2011; Prufer et al., 2015). Our survey involved a LiDAR technology called an Airborne Laser Swath Mapping (ALSM) system, consisting of a near infrared laser attached to a small airplane. The airplane's position and orientation in three-dimensional space was precisely measured with a combination of airborne and groundbased GPS and inertial navigation systems. The ALSM system captured a surface model of the study region that excluded ground cover once the data were appropriately filtered, leaving a model of the unobstructed surface. The preprocessing of the data and filtering of the surface model was done by the National Centre for Airborne Laser Mapping at the University of Houston (http://ncalm.cive.uh.edu). The surface model was then converted to a slope map and carefully surveyed by eye for Classic Maya settlements. Classic Maya settlements were identified and recorded using a modified version of the Xunantunich Archaeological Project classification system (Ashmore et al., 1994)—i.e., XAP typology (see Fig. 3). The XAP typology includes everything from isolated house-mounds (Type 1) to large patio groups (Type 5), and settlement units with structures exceeding 5 m in height (Type 7). It is based on the number of structures, their degree of formal arrangement (orthogonal or haphazard), and the height of the tallest structure in the cluster. However, the typology was not designed to account for larger centres—e.g. Caracol, a massive urban centre, would be a Type 7, greatly understating its size and complexity—so we provisionally added Types 8 and 9 for medium and large civic-ceremonial centres. Thus, sites in our survey were defined as any architectural feature, or features, that correspond to the XAP typology, meaning that some sites were individual domestic structures, others were clusters of domestic structures, and still others were civicceremonial centres. During the survey, we also gave every feature a confidence value (CV) to express our confidence about the visual identification of archaeological features. CV-1 indicates no confidence about the anthropogenic nature of the feature. CV-2 indicates confidence that the feature is anthropogenic, but no confidence about the given XAP classification. And CV-3 indicates confidence about both the anthropogenic nature of the feature and its XAP classification. The LiDAR analysis began with a blind test to estimate the accuracy of conducting a survey using LiDAR data. Dan Savage, who at the time was unfamiliar with the study region, searched the LiDAR surface for anthropogenic features within the study area. The full results of his analysis are currently being written up. Fig. 4 is a good example of a clear site from the LiDAR data. The sites and their CVs were recorded in a point data Shapefile using GRASS GIS. Savage then systematically compared his LiDAR survey results to the results of previous pedestrian surveys in three areas. The pedestrian survey areas encompassed the Classic Maya centres of Minanha and Waybil, and a nearby region called the Contreras Valley Zone (Lamoureux-St-Hilaire et al., 2015; McCane et al., 2010). After comparing his LiDAR survey results to the field reports from the ground surveys, he determined that he had approximately an 80% accuracy rate for large sites—i.e., XAP Type V and higher. His success rate for smaller sites was much lower, around 16%. The larger sites consisted of taller structures arranged orthogonally, making them much more clearly visible since they stood out as anthropogenic. This relationship between success rates and XAP types was also reflected in CVs, with higher confidence assigned to larger sites. Thus, the blind test indicated that the LiDAR survey results are skewed in favour of larger sites. However, since the LAMAP prediction we were testing was developed using large Classic Maya civic-ceremonial centres, we reasoned that our analyses would be accurate despite the bias. In total, Savage identified 8513 sites using the LiDAR images covering approximately 400 km2 of previously unsurveyed terrain. Of these, 4930 were given CVs of 2 or 3 (see Fig. 5), meaning he was confident that the features were anthropogenic and fairly confident in his assessment of their XAP types. XAP types I–III account for 1732 of the CV-2 and 64 W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 I II III IV V VI VII Isolated mound (< 2 m high) 2–4 mounds arranged informally (all < 2 m high) 2–4 mounds orthogonally arranged (< 2 m high) 5 or more mounds arranged informally (< 2 m high) 5 or more mounds with at least 2 arranged orthogonally (< 2 m high) 1 or more mounds arranged in any way (at least one > 2 m high) 1 or more mounds arranged in any way (at least one > 5 m high) XAP Typology. Fig. 3. XAP typology. CV-3 sites, and XAP types V–VII account for the remaining 314 sites. Of the sites that fell within the previously surveyed areas around Minanha, Waybil, and the Contreras Valley, 45 were previously unrecorded. They consisted mainly of XAP types I and II with 7 formally arranged type III groups. Unfortunately, given concerns about looting, a ongoing problem for heritage management in Belize, we cannot publish the detailed survey data, and interested parties should contact the authors directly. Once the visual survey was complete, we wanted to know what some of the unrecorded LiDAR features were like on the ground so that we could further assess the meaning of the CVs. To make this comparison, we used Barry and Cheong's (2012) field observations. Some of their field observations included areas Savage surveyed that were not also covered by the three other pedestrian surveys. Thus, Barry and Cheong's (2012) observations complemented our confidence assessment of the LiDAR survey data because we could compare their field notes and photos to the LiDAR results, at least subjectively. The field team's observations supported Savage's finding that the sites assigned CV-1—i.e., the lowest confidence level—were questionably anthropogenic while the sites assigned CV-2 were more likely archaeological and those assigned to CV-3 appeared to be accurately identified as both archaeological and from the Classic Maya period, corresponding to their XAP classification. Consequently, we excluded the CV-1 sites from further analysis. In the second stage of our analysis, we compared the survey point data to the LAMAP model from Carleton et al. (2012). This stage proceeded in several steps and involved a combination of GRASS GIS and the R statistical software environment. In the first step, we binned the LAMAP data into 5 archaeological potential classes using quantiles of the LAMAP values. The quantiles were computed in GRASS GIS and the bin edges can be seen in Table 1. Then, each archaeological feature identified in the LiDAR survey was assigned a LAMAP class based on its location in the LAMAP model. In the next step, we tested for differences in LAMAP classes between CV-2 and 3 sites. To do this, we used the Mann-Whitney U test. Since the LAMAP values are spatially autocorrelated, we generated the null distribution for the test with a bootstrap (Dodge, 2008; Efron, 1979). The null distribution is the expected density of the Mann-U statistic if there was no difference between the LAMAP values at CV-2 and CV-3 sites. If there was a difference, then the Mann-U statistic calculated using the empirical samples of CV-2 and CV-3 sites would be outside the 95% confidence region of the null distribution. To generate the distribution, we performed a bootstrap by pooling the CV-2 and CV-3 sites together and randomly sampling without replacement from the pooled sites. Each iteration of the bootstrap involved drawing two samples with sizes Fig. 4. Example of a site clearly visible on the LiDAR surface. Fig. 5. Map showing LAMAP model with CV2 and CV3 sites. W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 equal to the number of CV-2 and CV-3 sites respectively and then calculating the Mann-Whitney U statistic with R's ‘wilcox.test’ function. The bootstrap simulation was run for 10,000 iterations. It allowed us to overcome the problem with spatial autocorrelation that might have invalidated the test. Our results showed that there was no difference in LAMAP classes between the CV-2 and CV-3 sites, meaning that we could use either group to evaluate the LAMAP model (p = 0.4; see Fig. 6). So, given the additional confidence associated with the CV-3 sites and the increased likelihood that they were comparable to the sites used to build the LAMAP prediction, we focused our attention on CV-3 sites for the rest of the analysis. In the third step, we tested whether the CV-3 sites identified in the LiDAR survey could be distinguished from random locations in the LAMAP model. We reasoned that if they could not be distinguished, then the LAMAP model was performing no better than random chance—i.e., the distribution of identified CV-3 sites was random with respect to the LAMAP prediction. To conduct the test, we used the same bootstrap Mann-Whitney U test as before. We randomly selected 2608 point-locations on the LAMAP surface, corresponding to the number of CV-3 sites. We then extracted the LAMAP classes from these point locations and pooled these data with the LAMAP values from the CV-3 sites and ran the Mann-Whitney U bootstrap procedure as before. Our results showed that there was a statistically significant difference between the surveyed locations and the random sample, indicating that the LAMAP model was better on average than random chance (p ≤ 0.01). We did not provide a supporting image like the one in Fig. 6 because the value of the test statistic could not be plotted with the null distribution—it would be well off the page to right. In the penultimate step, we borrowed a concept from Economics called stochastic dominance to determine whether the LAMAP prediction consistently outperformed random chance. Stochastic dominance is commonly used when comparing lotteries to determine whether one gamble is better than another given the distribution of probable outcomes for each lottery (Hadar and Russell, 1969). If one lottery stochastically dominates another, then the dominant one should always be Fig. 6. Mann-Whitney U Bootstrap Test showing no statistically significant difference between the LAMAP values for CV-2 and CV-3 sites. The grey region is the bootstrapped distribution of the Mann-U statistic and the vertical line shows the value of the Mann-U test statistic for the empirical CV-2 and CV-3 site data, indicating that the CV-2 and CV-3 sites are statistically indistinguishable in terms of LAMAP class (i.e., quantile). We added the p-value for the test to the plot as well. The light grey area shows the 95% CI for the bootstrapped test statistic distribution. 65 chosen to maximize returns. To cast the concept in terms of archaeological potential, consider a situation in which you have to decide where to invest resources for finding unrecorded archaeological sites. Given two options, either using the LAMAP model or randomly selecting locations on the landscape, you would want to know which option was going to return the most sites. To evaluate the LAMAP performance, we used two lotteries: one in which we randomly sample the landscape and another in which we only sample from CV-3 site locations, tallying the number of times each LAMAP class is drawn in each lottery. To express the probabilities for all outcomes of the two lotteries, we calculated the empirical cumulative distribution functions (CDFs) for the CV-3 sites and randomly selected locations (Fig. 7). The CDFs show the probability of obtaining a LAMAP class less than or equal to each value shown on the x-axis. If the LAMAP prediction performed no better than random chance, we would expect the CV-3 sites to be distributed randomly with respect to LAMAP classes and the CDFs would be equal. On the other hand, if drawing from the pool of CV-3 sites was more likely to yield higher LAMAP classes we would expect the CV-3 distribution to stochastically dominate the random CDF, indicating that greater numbers of CV-3 sites are associated with higher LAMAP values. In terms of the plot depicted in Fig. 7, stochastic dominance would imply that the CV-3 CDF should be lower than or equal to the random CDF. As Fig. 7 shows, the CV-3 CDF stochastically dominates the random CDF. In fact, the comparison indicates first order stochastic dominance, an especially strong type of dominance in which the CV-3 CDF is always more likely to yield higher LAMAP values than randomly selected locations. Thus, the LAMAP model appears to always a better bet for finding more sites. In the last step, we sought to determine the relationship between LAMAP potential values and the numbers of sites actually identified by our survey. Together, the two tests of randomness could only tell us that the LAMAP model consistently outperformed random chance, not how well the model performed in easily interpretable terms. A more important test would be whether real, substantial gains in the numbers of sites found would come from using the model. To perform the test, we estimated a linear regression in R between the number of sites identified in the survey and the corresponding LAMAP classes associated Fig. 7. Showing the first-order stochastic dominance of the CV-3 sites over random locations. The light grey area represents the cumulative distribution function (CDF) of the CV-3 LAMAP classes and the dark grey area shows the CDF of the LAMAP classes of randomly selected points. 66 W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 with each site. Our results showed that the LAMAP model performed very well. There is a clear relationship between LAMAP classes and the quantity of sites identified by our survey (see Fig. 8). The correlation is very strong (R2 = 0.96; p ≤ 0.01), and the linear model indicates that each increase in LAMAP class was related to an additional 138 sites, on average. Between the lowest LAMAP class and the highest, there is an increase of approximately 550 sites. 4. Discussion and conclusions Several predictive tools for estimating archaeological potential are available and have been used in academic research with varying degrees of success. However, our recently developed technique called the LAMAP method has several advantages over the other methods. At the time of its inception, we used a case study involving Classic Maya sites to demonstrate the technique. In that study we created a LAMAP prediction for a large portion of the Vaca Plateau around the Classic Maya centre of Minanha, located in west-central Belize. The assessment was validated by withholding a small subset of known sites before calculating the model and then using that subset to test the model's predictions. While this initial test suggested that the LAMAP method was performing well, the test was not comprehensive. So, we conducted a more comprehensive test by surveying the Minanha study area with a combination of remotely sensed LiDAR data and fieldwork. By comparing the newly acquired survey data to the LAMAP prediction, we tested the hypothesis that the LAMAP method provides accurate predictions of archaeological potential. The hypothesis contends that there should be a statistically significant relationship between the LAMAP potential assessment and the number of sites found. Our results indicate that a strong positive relationship exists. The relationship suggests an archaeologist looking to find sites in the study area would be able to find as many as 550 additional sites by concentrating on areas predicted by the LAMAP model to have high archaeological potential. These findings are statistically significant compared to randomness. Thus, the hypothesis seems to hold after a comprehensive test. Fig. 8. Showing a histogram of the frequencies of CV-3 sites in each LAMAP class. Overlaid is a scatter plot and the linear model for the regression of CV-3 sites on LAMAP class. The regression shows a positive linear relationship between the quantity of CV-3 sites and the LAMAP, indicating that the LAMAP model accurately reflects archaeological potential in the study area. Our findings have three main implications. One is that the LAMAP method can be used to direct survey efforts in rugged, difficult to access terrain with dense forest cover—conditions that often make extensive surveys prohibitively expensive. The method is a cost-effective, accurate way of assessing archaeological potential. Like other predictive methods, LAMAP predictions use remotely sensed landscape and environmental data. With many datasets now freely available online, including high-resolution multi-spectral data from NASA satellite platforms like ASTER, the cost of calculating a LAMAP prediction is miniscule compared to the costs of extensive pedestrian surveys. And although the model does not replace fieldwork, it does allow archaeologists to use their resources more effectively, improving the return on research investment. Another implication is that the variability in land-use of the North Vaca Plateau by the Classic Maya might be captured surprisingly well by a small number of landscape variables. The LAMAP model of the region surrounding Minanha involved only four landscape variables, namely slope, ruggedness, distance to water, and soil arability. Even with this small number of variables, the model produced potential estimates that strongly correlate with archaeological resource quantities. This correlation suggests that the landscape variables involved in the prediction were very important for the Classic Maya. We suspect that an improved soil arability dataset would enhance the predictions further because of the Classic Maya's dependence on agriculture. We also think that future research should look into including other types of variables like visibility models (e.g., Barry, 2014; Llobera, 2003) or sociocultural parameters such as the presence of ritually important features like caves (Moyes et al., 2009). Future efforts should also involve comparing different sets of variables to find the optimal combinations. Such research would allow us to test the individual importance of different landscape variables to the Classic Maya and shed additional light on past human–landscape interaction. A third implication is that the theory underlying the LAMAP method might be an accurate representation of human land-use decisions. To reiterate, the LAMAP theory states that when people plan to engage in a given activity on the landscape, they decide where to do it by comparing candidate locations for the activity to one or more sociocultural archetypes. The archetypes describe the landscape characteristics appropriate for that activity—e.g., the best place to camp, or found a settlement. They are the product of long-term human–landscape interaction and they are likely passed on through social learning and maintained in cultural memory. For example, people will tend to camp where others in their communities camp or have camped before because the landscape characteristics associated with others' campsites have turned out to be optimal for camping. Given the success of the LAMAP method so far, it seems that this theory might go some way toward explaining how humans actually made land-use decisions in the past. Still, many other tests are required, such as determining whether there is evidence for the inheritance of archetypes or testing the theory using contemporary people, indicating profitable avenues for future research. The last implication is that the LAMAP technique appears to be useful for managing heritage resources. Archaeological resources are under pervasive threat all over the world. Climate change, modern industry, and tourism disturb many sites and their surrounding landscapes, threatening their integrity and precipitating a loss of knowledge about the past (e.g., Durand et al., 2010; Howard et al., 2008; Hutchinson, 1996; Murphy et al., 2009). And while the loss of knowledge is troubling, it is often difficult or impossible to avoid because the processes that lead to it are necessary or unavoidable. Necessary damages result from economic activity, including natural resource extraction and infrastructure development, while unavoidable damages result from natural degradation and climate change, involving processes like sea level rise, extreme flooding, or erosion. As a result, the degradation of archaeological resources is inevitable. Unfortunately, managing the damage is challenging because the finances, time, and expertise required are not W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 always available for preventing the damage or recording sites before they are lost, problems that are compounded for unrecorded archaeological resources because we cannot manage resources when we are not aware of them. Ultimately, the problem is that the losses can easily outstrip our capacity to mitigate them. With the integrity of the archaeological record under threat globally from human and natural degradation, the potential loss of information about the past is worrying. And surveying previously unexplored regions to determine which areas are under the greatest threat of degradation is often prohibitively expensive and time consuming, as our field survey efforts demonstrated. Thus, it is important to develop predictive tools that can help to manage archaeological resources. The LAMAP method appears to have substantial potential in that regard. Using its predictions, archaeologists can help governments to better manage heritage resources, especially in places where resources are limited and the potential loses are significant. Importantly, the LAMAP method can also help us understand past human land-use behaviour by reconstructing the decision-making processes of past people, moving us beyond prediction and toward understanding. Acknowledgements The authors express their gratitude to the Institute of Archaeology in Belize, and specially Drs. Jaime Awe and John Morris, for their staunch support of our various research endeavours. The LiDAR study was graciously funded through a consortium grant provided by the Alphawood Foundation. Much of the reconnaissance, mapping, and excavations at Minanha was conducted under a series of grants from the Social Sciences and Humanities Research Council of Canada. Carleton also received financial support from the Human Evolutionary Studies Program and the Department of Archaeology at Simon Fraser University and would like to thank his supervisor, Prof. Mark Collard, for his comments on an earlier draft of this paper. Finally, we give thanks to all of the members of the Social Archaeology Research Program who’s hard work have brought the past of the North Vaca Plateau to light. References Alexakis, D., Sarris, A., Astaras, T., Albanakis, K., 2011. Integrated GIS, remote sensing and geomorphologic approaches for the reconstruction of the landscape habitation of Thessaly during the Neolithic period. J. Archaeol. Sci. 38, 89–100. Andrica, D., Andreescu, T., 2009. Number Theory. Birkhäuser Boston, Boston. Ashmore, W., Connell, S., Ehret, J., Gifford, C., Neff, T., Vandenbosch, J., 1994. The Xunantunich settlement survey. In: Leventhal, R. (Ed.), Xunantunich Archaeological Project: 1994 Field Season. Xunantunich Archaeological Project, pp. 248–288. Barbet-Massin, M., Jiguet, F., Albert, C.H., Thuiller, W., 2012. Selecting pseudo-absences for species distribution models: how, where and how many? Methods Ecol. Evol. 3: 327–338. http://dx.doi.org/10.1111/j.2041-210X.2011.00172.x. Barry, J., 2014. A GIS-based spatial analysis of visibility and movement using the ancient Maya center of Minanha, Belize. (Unpublished MA Thesis). Trent University, Peterborough, ON. Barry, J., Cheong, K., 2012. Over the hills and through the bush: reconnaissance in the Minanha hinterlands during the 2012 field season. In: Iannone, G., Awe, J., Schwake, S., Hills, K. (Eds.), Archaeological Investigations in the North Vaca Plateau, Belize: Progress Report of the Fourteenth (2012) Field Season. Department of Anthropology, Trent University, Peterborough, ON. Bevan, A., Conolly, J., 2002. GIS, archaeological survey, and landscape archaeology on the island of Kythera, Greece. J. Field Archaeol. 29:123–138. http://dx.doi.org/10.2307/ 3181488. Bewley, R.H., Crutchley, S.P., Shell, C.A., 2005. New light on an ancient landscape: lidar survey in the Stonehenge World Heritage Site. Antiquity 79, 636–647. Binford, L.R., 1980. Willow smoke and dogs' tails: hunter-gatherer settlement systems and archaeological site formation. Am. Antiq. 4–20. Binford, L.R., 1982. The archaeology of place. J. Anthropol. Archaeol. 1, 5–31. Bird, D.W., O'Connell, J.F., 2006. Behavioral ecology and archaeology. J. Archaeol. Res. 14: 143–188. http://dx.doi.org/10.1007/s10814-006-9003-6. Carleton, W.C., Conolly, J., Ianonne, G., 2012. A locally-adaptive model of archaeological potential (LAMAP). J. Archaeol. Sci. 39:3371–3385. http://dx.doi.org/10.1016/j.jas.2012. 05.022. Challis, K., 2006. Airborne laser altimetry in alluviated landscapes. Archaeol. Prospect. 13: 103–127. http://dx.doi.org/10.1002/arp.272. Chase, A.F., Chase, D.Z., Weishampel, J.F., Drake, J.B., Shrestha, R.L., Slatton, K.C., Awe, J.J., Carter, W.E., 2011. Airborne LiDAR, archaeology, and the ancient Maya landscape at Caracol, Belize. J. Archaeol. Sci. 38:387–398. http://dx.doi.org/10.1016/j.jas.2010.09.018. 67 Chase, A.F., Chase, D.Z., Fisher, C.T., Leisz, S.J., Weishampel, J.F., 2012. Geospatial revolution and remote sensing LiDAR in Mesoamerican archaeology. Proc. Natl. Acad. Sci. U. S. A. 109:12916–12921. http://dx.doi.org/10.1073/pnas.1205198109. Chefaoui, R.M., Lobo, J.M., 2008. Assessing the effects of pseudo-absences on predictive distribution model performance. Ecol. Model. 210:478–486. http://dx.doi.org/10. 1016/j.ecolmodel.2007.08.010. Conolly, J., Lake, M., 2006. Geographical Information Systems in Archaeology. Cambridge University Press, Cambridge. Conolly, J., Manning, K., Colledge, S., Dobney, K., Shennan, S., 2012. Species distribution modelling of ancient cattle from early Neolithic sites in SW Asia and Europe. The Holocene 22:997–1010. http://dx.doi.org/10.1177/0959683612437871. Demján, P., Dreslerová, D., 2016. Modelling distribution of archaeological settlement evidence based on heterogeneous spatial and temporal data. J. Archaeol. Sci. 69: 100–109. http://dx.doi.org/10.1016/j.jas.2016.04.003. Dodge, Y., 2008. The Concise Encyclopedia of Statistics. Springer. Dorshow, W.B., 2012. Modeling agricultural potential in Chaco Canyon during the Bonito phase: a predictive geospatial approach. J. Archaeol. Sci. 39:2098–2115. http://dx.doi. org/10.1016/j.jas.2012.02.004. Durand, J.F., Meeuvis, J., Fourie, M., 2010. The Threat of Mine Effluent to the UNESCO Status of the Cradle of Humankind World Heritage Site. Durand-Dastès, F., Favory, F., Fiches, J., Mathian, H., Pumain, D., Raynaud, C., Sanders, L., Van der Leeuw, S., 1998. Des oppida aux métropoles. Archéologues et géographes en vallée du Rhône. Anthropos. Ebert, D., 2004. Applications of archaeological GIS. Can. J. Archaeol. 28, 319–341. Efron, B., 1979. Bootstrap methods: another look at the jackknife. Ann. Stat. 7:1–26. http://dx.doi.org/10.1214/aos/1176344552. Fernandes, R., Geeven, G., Soetens, S., Klontza-Jaklova, V., 2011. Deletion/substitution/addition (DSA) model selection algorithm applied to the study of archaeological settlement patterning. J. Archaeol. Sci. 38:2293–2300. http://dx.doi.org/10.1016/j.jas.2011. 03.035. Graves, D., 2011. The use of predictive modelling to target Neolithic settlement and occupation activity in mainland Scotland. J. Archaeol. Sci. 38:633–656. http://dx.doi.org/ 10.1016/j.jas.2010.10.016. Hadar, J., Russell, W.R., 1969. Rules for ordering uncertain prospects. Am. Econ. Rev. 59, 25–34. Hastie, T., Fithian, W., 2013. Inference from presence-only data; the ongoing controversy. Ecography 36:864–867. http://dx.doi.org/10.1111/j.1600-0587.2013.00321.x. Hitchings, P., Bikoulis, P., Edwards, S., Banning, E.B., 2016. Predict and confirm: Bayesian survey and excavation at three candidate sites for late Neolithic occupation in Wadi Quseiba, Jordan. In: Campana, S., Scopigno, R., Carpentiero, G., Cirillo, M. (Eds.), Proceedings of the 43rd Annual Conference on Computer Applications and Quantitative Methods in Archaeology. Presented at the Computer Applications and Quantitative Methods in Archaeology. Archaeopress, Toronto, pp. 605–612. Howard, A.J., Challis, K., Holden, J., Kincey, M., Passmore, D.G., 2008. The impact of climate change on archaeological resources in Britain: a catchment scale assessment. Clim. Chang. 91:405–422. http://dx.doi.org/10.1007/s10584-008-9426-9. Howey, M.C.L., 2015. Geospatial landscape permeability modeling for archaeology: a case study of food storage in northern Michigan. J. Archaeol. Sci. 64:88–99. http://dx.doi. org/10.1016/j.jas.2015.10.007. Hutchinson, G., 1996. Threats to underwater cultural heritage. The problems of unprotected archaeological and historic sites, wrecks and objects found at ea. Marine Policy, the Protection of the Underwater Cultural Heritage Vol. 20:pp. 287–290. http://dx.doi. org/10.1016/0308-597X(96)00017-6. Iannone, G., 2005. The rise and fall of an ancient Maya petty royal court. Lat. Am. Antiq. 16:26–44. http://dx.doi.org/10.2307/30042485. Iannone, G., 2010. Collective memory in the Frontiers: a case study from the ancient Maya Centre of Minanha, Belize. Anc. Mesoam. 21:353–371. http://dx.doi.org/10.1017/ S09565361100002508. Iannone, G., Chase, A., Chase, D., Awe, J., Moyes, H., Brook, G., Polk, J., Webster, J., Conolly, J., 2014a. An archaeological consideration of long-term socioecological dynamics on the Vaca Plateau, Belize. The Great Maya Droughts in Cultural Context: Case Studies in Resilience and Vulnerability. University Press Of Colorado. Iannone, G., Prufer, K., Chase, D.Z., 2014b. Resilience and vulnerability in the Maya hinterlands. Archeol. Pap. Am. Anthropol. Assoc. 24:155–170. http://dx.doi.org/10.1111/ apaa.12035. Kanter, J., 2008. The archaeology of regions: from discrete analytical toolkit to ubiquitous spatial perspective. J. Archaeol. Res. 16, 37–81. Kvamme, K., 1988. Development and testing of quantitative models. In: Judge, W., Sebastian, L. (Eds.), Quantifying the Present and Predicting the Past: Theory, Method and Application of Archaeological Predictive Modeling. US Bureau of Land Management, Department of the Interior, US Government Printing Office, Washington, DC. Lamoureux-St-Hilaire, M., Macrae, S., McCane, C., Parker, E., Iannone, G., 2015. The last groups standing: living abandonment at the ancient Maya Center of Minanha, Belize. Lat. Am. Antiq. 26, 550–569. Llobera, M., 2003. Extending GIS-based visual analysis: the concept of visualscapes. Int. J. Geogr. Inf. Sci. 17:25–48. http://dx.doi.org/10.1080/713811741. McCane, C., Macrae, S., Iannone, G., 2010. A consideration of the spatial arrangement of settlement groups and terraces in Contreras, Minanha, Belize. Res. Rep. Belizean Archaeol. 7, 141–151. McCoy, M.D., Ladefoged, T.N., 2009. New developments in the use of spatial technology in archaeology. J. Archaeol. Res. 17:263–295. http://dx.doi.org/10.1007/s10814-0099030-1. Mehrer, M.W., Wescott, K.L., 2006. GIS and Archaeological Site Location Modeling. CRC Press. Moyes, H., Awe, J.J., Brook, G.A., Webster, J.W., 2009. The ancient Maya drought cult: late classic cave use in Belize. Lat. Am. Antiq. 20, 175–206. 68 W.C. Carleton et al. / Journal of Archaeological Science: Reports 11 (2017) 59–68 Murphy, P., Thackray, D., Wilson, E., 2009. Coastal heritage and climate change in England: assessing threats and priorities. Conserv. Manag. Archaeol. Sites 11:9–15. http://dx.doi.org/10.1179/135050309X12508566208281. Oonk, S., Spijker, J., 2015. A supervised machine-learning approach towards geochemical predictive modelling in archaeology. J. Archaeol. Sci. 59:80–88. http://dx.doi.org/10. 1016/j.jas.2015.04.002. Phillips, S.J., Anderson, R.P., Schapire, R.E., 2006. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190:231–259. http://dx.doi.org/10.1016/j. ecolmodel.2005.03.026. Prufer, K.M., Thompson, A.E., Kennett, D.J., 2015. Evaluating airborne LiDAR for detecting settlements and modified landscapes in disturbed tropical environments at Uxbenká, Belize. J. Archaeol. Sci. 57:1–13. http://dx.doi.org/10.1016/j.jas.2015.02.013. R Core Team, 2016. R: A Language and Environment for Statistical Computing. Rua, H., 2009. Geographic information systems in archaeological analysis: a predictive model in the detection of rural roman Villae. J. Archaeol. Sci. 36, 224–235. Schwake, S.A., Iannone, G., 2010. Ritual remains and collective memory: Maya examples from west Central Belize. Anc. Mesoam. 21:331–339. http://dx.doi.org/10.1017/ S0956536110000283. Smith, E.A., Winterhalder, B., 1992. Evolutionary Ecology and Human Behavior. Transaction Publishers. Swanson, S., 2003. Documenting prehistoric communication networks: a case study in the Paquimé polity. Am. Antiq. 68, 753–767. Tallavaara, M., Seppä, H., 2012. Did the mid-Holocene environmental changes cause the boom and bust of hunter-gatherer population size in eastern Fennoscandia? The Holocene 22:215–225. http://dx.doi.org/10.1177/0959683611414937. Archéologie et systèmes socio-environnementaux. In: Van der Leeuw, S., Favory, F., Fiches, J. (Eds.), Études multiscalaires sur la vallée du Rhône dans le programme Archaeomedes. CNRS Éditions, Paris. Verhagen, P., 2007. Case Studies in Archaeological Predictive Modelling. Amsterdam University Press. Verhagen, P., Whitley, T.G., 2011. Integrating archaeological theory and predictive modeling: a live report from the scene. J. Archaeol. Method Theory 19:49–100. http://dx. doi.org/10.1007/s10816-011-9102-7. Verhagen, P., Nuninger, L., Tourneux, F., Bertoncello, F., Jeneson, K., 2013. Introducing the human factor in predictive modelling: a work in progress. In: Earl, G., Sly, T., Chrysanthi, A., Murrieta-Flores, P., Papadopoulos, C., Romanowska, C., Wheatley, D. (Eds.), Archaeology in the Digital Era. Papers from the 40th Annual Conference of Computer Applications and Quantitative Methods in Archaeology (CAA). Amsterdam University Press, pp. 379–388. Wheatley, D., Gillings, M., 2002. Spatial Technology and Archaeology. Taylor and Francis, London. William, B., D'Errico, F., Dibble, H.L., Krishtalka, L., West, D., Olszewski, D., Peterson, A.T., David, G.A., Gillam, G.C., Montet-White, A., Crucifix, M., Marean, C.W., SánchezGoñi, M.-F., Wohlfarthh, B., Vanhaeren, M., 2006. Eco-cultural niche modeling: new tools for reconstructing the geography and ecology of past human populations. Palaeontol. Anthropol. 68–83.